DESIGN  OF  A  LIGHTER  THAN  AIR  VEHICLE  THAT  ACHIEVES  POSITIVE 
BUOYANCY  IN  AIR  USING  A  VACUUM 


THESIS 


Trent  T.  Metlen,  Captain,  USMC 


AF  IT -EN  Y  -13-J-02 


DEPARTMENT  OF  THE  AIR  FORCE 
AIR  UNIVERSITY 

AIR  FORCE  INSTITUTE  OF  TECHNOLOGY 


Wright-Patterson  Air  Force  Base,  Ohio 


APPROVED  FOR  PUBLIC  RELEASE;  DISTRIBUTION  UNLIMITED 


The  views  expressed  in  this  thesis  are  those  of  the  author  and  do  not  reflect  the  official 
policy  or  position  of  the  United  States  Air  Force,  Department  of  Defense,  or  the  United 
States  Government. 


AF  IT -EN  Y  -13-J-02 


DESIGN  OF  A  LIGHTER  THAN  AIR  VEHICLE  THAT  ACHIEVES  POSITIVE 
BUOYANCY  IN  AIR  USING  A  VACUUM 


THESIS 


Presented  to  the  Faculty 
Department  of  Aeronautics  and  Astronautics 
Graduate  School  of  Engineering  and  Management 
Air  Force  Institute  of  Technology 
Air  University 

Air  Education  and  Training  Command 
In  Partial  Fulfillment  of  the  Requirements  for  the 
Degree  of  Master  of  Science  in  Aeronautical  Engineering 


Trent  T.  Metlen,  BS 
Captain,  USMC 
June  2012 


APPROVED  FOR  PUBLIC  RELEASE;  DISTRIBUTION  UNLIMITED. 


AF  IT -EN  Y  -13-J-02 


DESIGN  OF  A  LIGHTER  THAN  AIR  VEHICLE  THAT  ACHIEVES  POSITIVE 
BUOYANCY  IN  AIR  USING  A  VACUUM 


Trent  T.  Metlen,  BS 
Captain,  USMC 


Approved: 


Dr.  Anthony  Palazotto  (Chainnan) 


Date 


Dr.  Jeremy  Agte,  Lt  Col,  USAF  (Member) 


Date 


Dr.  Mark  Reeder  (Member) 


Date 


Abstract 


Three  designs  for  a  Lighter  Than  Air  (LTA)  structure  that  achieve  positive 
buoyancy  using  a  vacuum  in  place  of  a  lifting  gas  were  proposed  and  evaluated.  Before 
the  first  human  flight  it  was  predicted  that  LTA  flight  was  possible  through  the  use  of  a 
light  weight  structure  maintaining  an  internal  vacuum.  Since  that  time  LTA  flight  has 
been  used  to  accomplish  various  missions  through  the  use  of  lifting  gases.  This  study  was 
conducted  in  response  to  an  anticipated  shortage  of  helium,  the  danger  of  hydrogen,  and 
the  possibility  of  using  LTA  vehicles  as  a  means  of  passenger  or  cargo  transportation. 

The  improvements  required  of  traditional  LTA  vehicles  to  become  viable  methods  of 
transportation  were  examined  and  found  to  be  buoyancy  control,  ground  handling,  and 
efficiency.  Efficiency  was  evaluated  using  von  Karman  efficiency  which  enabled 
comparison  between  vehicles  of  different  types  and  missions  based  on  velocity,  payload 
and  power  required.  Von  Karman  efficiency  was  also  used  to  compare  theoretical  LTA 
vehicles  of  constant  size  over  a  range  of  engine  sizes.  Vacuum  LTA  structures  were 
evaluated  based  on  their  ratio  of  structure  weight  to  the  weight  of  displaced  air  and 
termed  Weight/Buoyant  Force  (W/B)  where  a  W/B<1  corresponds  to  positive  buoyancy. 

Two  methods  of  geometrically  stiffening  a  sphere  were  investigated.  A  thin 
shelled  sphere  stiffened  with  an  isogrid  of  blade  type  stiffeners  using  Ultra  High  Modulus 
(UHM)  carbon  epoxy  was  predicted  to  give  a  W/B=  0.81  which  would  enable  LTA  flight. 
The  same  design  with  a  beryllium  skin  and  UHM  carbon  epoxy  stiffeners  was  predicted 
to  have  a  W/B=0.79.  A  geodesic  sphere  composed  of  a  frame  of  cylindrical,  UHM, 
carbon  epoxy,  pultruded  rods  with  a  thin  membrane  stretched  over  it  was  evaluated  using 


IV 


Finite  Element  Analysis  with  the  frame  modeled  as  clamped,  cylindrical  beam  elements. 
A  W/B=0.57  was  calculated  for  the  frame.  The  skin  was  predicted  to  add  approximately 
0.37  to  the  W/B  if  Zylon  was  used  to  reinforce  a  Mylar  membrane,  which  would  result  in 
a  structure  with  an  overall  W/B=0.94.  If  a  material  as  strong  as  graphene  were  used  for 
the  skin  the  skin  would  have  a  W/B=0.001,  resulting  in  an  overall  W/B=0.57.  These 
structures  were  optimized  using  a  non-linear  programming  optimization  routine  in 
Matlab. 

A  unique  LTA  vehicle  concept  composed  of  twin  counter-rotating  cylinders  that 
made  use  of  a  virtual  elliptical  body  to  reduce  drag  and  maintain  structural  integrity  with 
a  vacuum  was  also  investigated  and  found  to  be  infeasible. 
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DESIGN  OF  A  LIGHTER  THAN  AIR  VEHICLE  THAT  ACHIEVES  POSITIVE 
BUOYANCY  IN  AIR  USING  A  VACUUM 


I.  Introduction 


1.1  Motivation 

Lighter  Than  Air  (LTA)  vehicles  have  historically  used  lifting  gases  to  achieve 
positive  buoyancy  in  air.  Buoyancy  is  defined  as  the  force  exerted  on  an  object  by  the 
fluid  in  which  it  is  submerged.  This  force  is  due  to  the  difference  in  pressure  at  the  top  of 
an  object  vs.  the  bottom  of  the  object.  There  is  a  difference  in  pressure  because  of  the 
weight  of  the  column  of  fluid  surrounding  the  object.  Positive  buoyancy  is  defined  as  the 
weight  of  an  object  being  less  than  the  weight  of  the  fluid  it  displaces.  Negative  buoyancy 
is  defined  as  the  weight  of  an  object  being  more  than  the  weight  of  the  fluid  it  displaces, 
in  other  words,  it  tends  to  sink.  Neutrally  buoyant  means  the  object  is  the  same  weight  as 
the  fluid  it  displaces,  so  that  the  force  of  buoyancy  is  equal  to  the  object’s  weight  [1].  A 
vehicle  could  theoretically  achieve  positive  buoyancy  through  the  use  of  a  vacuum  in 
place  of  a  lifting  gas.  The  difficulty  of  designing  a  vacuum  LTA  structure  has  prevented 
it  from  being  used  in  the  past.  The  primary  advantage  of  a  vacuum  LTA  structure  would 
be  the  elimination  of  the  requirement  to  purchase  or  transport  a  lifting  gas. 

If  a  vacuum  could  be  used  to  replace  lifting  gas,  it  would  require  a  suitable 
mission  to  make  it  more  than  just  an  expensive  parlor  trick.  There  are  various  niche 
missions  currently  filled  by  LTA  vehicles,  but  they  are  currently  unable  to  compete  with 
the  mainstream  modes  of  transportation.  The  vacuum  LTA  vehicle  was  evaluated  to 
determine  missions  for  which  it  was  suited.  One  such  mission  motivated  the  first  sponsor 


1 


of  this  research  to  investigate  the  following  line  of  reasoning:  Small,  regional  airports  are 
largely  underused  while  major  airports  are  often  congested.  A  vehicle  capable  of 
efficiently  carrying  cargo  while  operating  out  of  small  airports  would  fill  a  niche 
transportation  market.  This  would  require  a  vehicle  that  could  land  on  short  runways  with 
large  payloads.  Could  the  vacuum  LTA  vehicle  possibly  be  a  solution  to  this  challenge? 

1.2  Background 

In  1670  Francesco  Lana  de  Terzi  published  a  book  titled,  “Prodromo  dell’ Arte 
Maestra”  in  which  he  proposed  the  design  of  a  LTA  vehicle  that  he  theorized  would 
attain  positive  buoyancy  in  air  using  a  vacuum  [2],  His  structure  resembled  a  ship  and 
was  suspended  from  four  thin-shelled  copper  spheres  as  shown  in  Figure  1  [3]  . 
Unfortunately  copper  does  not  have  a  high  enough  specific  stiffness  to  achieve  positive 
buoyancy  in  air  using  a  smooth  shell,  as  will  be  shown  in  the  theory  section  of  this  thesis. 
Nevertheless,  Lana  correctly  understood  that  the  principal  of  buoyancy  could  be  applied 
to  air  and  that  the  sphere  was  the  ideal  shape  for  a  LTA  vacuum  enclosure.  The  sphere  is 
the  ideal  shape  for  a  vacuum  enclosure  because  it  requires  the  minimum  thickness  to 
achieve  the  stiffness  required  to  prevent  buckling  of  any  shape.  It  also  has  the  minimum 
surface  area  per  volume  of  any  shape,  which  maximizes  buoyancy,  minimizes  pressure 
loading,  and  minimizes  weight. 
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Figure  1:  Lana’s  Flying  Boat  |3] 


During  the  course  of  this  research  it  was  not  uncommon  to  find  blogs  and  even 
published  papers  musing  about  the  possibility  of  vacuum  LTA  structures.  Most  of  the 
authors  did  not  apply  engineering  principles  to  the  problem.  The  most  common 
misconception  was  that  if  a  sphere  was  made  big  enough  the  problem  would  get  easier 
because  of  the  squared  cube  law.  The  idea  behind  the  squared  cube  law  is  that  the  area  of 
the  sphere  grows  at  a  rate  of  radius  squared,  while  the  displaced  volume  of  air  grows  at  a 
rate  of  radius  cubed.  This  is  true,  but  the  fallacy  of  this  argument  is  in  not  understanding 
that  the  thickness  of  the  sphere  must  grow  as  well.  The  tendency  for  a  thin  shelled  sphere 
to  buckle  under  external  pressure  is  a  function  of  the  thickness  of  the  shell  as  shown  by 
equation  (14)  [4,  p.  1292],  To  prevent  buckling,  the  thickness  must  grow  proportionally 
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with  the  surface  area  of  the  sphere,  so  that  no  matter  how  big  it  is,  the  weight  of  the 
sphere  remains  proportional  to  the  weight  of  the  air  displaced. 

Akhmeteli  and  Gavrilin  claimed  that  a  vacuum  LTA  structure  could  be 
constructed  using  layered  shell  sandwich  construction  [5].  Their  approach  to  the  vacuum 
sphere  formed  the  basis  for  equation  (23)  which  not  only  shows  why  copper  cannot  be 
used  to  fonn  a  thin  shelled,  positively  buoyant  vacuum  sphere,  but  why  no  material 
available  today  could  either.  Only  a  geometrically  stiffened  shell,  such  as  the  sandwich 
construction  investigated  by  Akhmeteli  and  Gavrilin,  could  possibly  achieve  this  feat. 

1.3  Research  Objective 

For  a  theoretical  vehicle  to  be  operationally  viable,  it  must  achieve  perfonnance 
comparable  to  or  better  than  existing  vehicles  with  similar  mission  types.  The  C-27J  was 
designed  to  transport  cargo  to  small  airports  relatively  efficiently.  One  purpose  of  this 
thesis  was  to  detennine  whether  a  vacuum  LTA  vehicle  could  effectively  perform  the 
same  mission  as  the  C-27J.  Therefore  it  became  necessary  to  evaluate  the  C-27J 
capabilities  and  compare  it  to  any  proposed  vacuum  LTA  vehicle  design.  The  C-27J 
specifications  are  listed  in  Table  1,  [6], 
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Table  1:  C-27J  Specifications  [7] 


Dimensions 


Span 

28.70  m 

94.16ft 

Length 

22.70  m 

74.48  ft 

Height 

9.65  m 

31.6  ft 

Wing  area 

81.94  m 

882  ft 

Weights 


Normal  Landing  Weight 

27,500  kg 

60,627  lb 

Maximum  Takeoff 

31,800  kg 

70,107  lb 

Power  Plant 


Engine 

Rolls-Royce  AE  2100-D2 

Engine  Power 

4637  Shaft  Horse  Power 

Propeller 

Dowty  R-391  six-blade 

Performance  (clean,  ISA) 

Maximum  Service  Ceiling 

9,144  m 

30,000ft 

Maximum  Cruise  Speed 

315  KTAS 

Landing  Ground  Roll 

340  m 

1,115  ft 

(Maximum  Landing  Weight,  Sea  Level) 


LTA  vehicles  are  generally  stable,  limited  in  speed,  have  long  loiter  times,  large 
payloads,  and  poor  handling  qualities  in  the  landing  environment.  These  characteristics 
have  largely  relegated  them  to  niche  markets  such  as  surveillance,  communication  relays, 
and  advertisement.  Table  2  lists  key  questions  that  must  be  answered  in  order  to 
determine  whether  a  vacuum  LTA  vehicle  is  a  feasible  concept. 
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Table  2:  Research  Questions 


1 

Can  positive  buoyancy  be  achieved  by  a  structure  in  air  using  a  vacuum? 

2 

Is  a  vacuum  LTA  structure  a  feasible  platform  for  any  existing  missions? 

3 

What  is  the  ideal  internal  pressure  of  a  partial  vacuum  LTA  structure? 

4 

What  is  the  ideal  operating  altitude  of  a  vacuum  LTA  structure? 

5 

Are  there  viable  geometric  shapes  for  a  vacuum  LTA  structure? 

6 

Do  construction  techniques  exist  which  can  construct  the  required 
geometric  shapes? 

7 

What  are  the  critical  material  properties  for  a  vacuum  LTA  structure? 

8 

Do  materials  exist  from  which  a  vacuum  LTA  structure  can  be 
constructed? 

1.4  Thesis  Overview 

Is  LTA  flight  possible  using  a  vacuum?  During  the  course  of  this  investigation 
von  Kannan’s  efficiency  was  used  to  compare  theoretical  vehicle  design  spaces  to 
existing  vehicles.  A  unique  vehicle  consisting  of  counter-rotating  cylinders  that  take 
advantage  of  a  virtual  elliptical  body  to  minimize  drag  was  proposed  and  evaluated  for  its 
energetic  feasibility.  A  program  was  developed  in  Matlab  that  optimized  the  geometry  of 
a  geodesic  sphere  through  the  use  of  built  in  optimization  functions  and  finite  element 
analysis.  All  of  these  investigations  attempted  to  answer  the  research  questions  posed  in 
Table  2. 

The  literature  review  section  examines  historical  LTA  vehicles  to  detennine  the 
criterion  for  a  successful  LTA  vehicle.  It  further  defines  the  problem  of  constructing  a 
vacuum  sphere  through  the  study  of  thin  shelled  spheres  under  pressure.  Finally,  it 
introduces  the  two  methods  of  geometrically  improving  upon  a  thin  shelled  sphere  that 
were  used  in  this  study.  These  are  stiffening  of  a  thin  shell  using  blade  type  stiffeners 
arranged  in  an  isogrid  and  a  geodesic  sphere  of  cylindrical  beams  supporting  a  thin 
membrane. 
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The  methodology  section  introduces  the  counter-rotating  cylinder  LTA  vehicle 
and  explains  the  methods  used  to  evaluate  its  feasibility.  It  explains  the  methods  used  to 
evaluate  an  isogrid  design.  Finally,  it  describes  how  the  geometry  of  a  geodesic  sphere 
was  calculated,  evaluated,  and  optimized. 

The  results  section  gives  an  answer  to  the  viability  of  each  of  the  proposed 
methods  for  achieving  positive  buoyancy  using  a  vacuum.  The  vehicle  design  spaces  are 
plotted  on  a  von  Kannan  efficiency  graph.  The  design  spaces  of  the  isogrid  stiffened 
sphere  and  geodesic  spheres  are  graphed  for  a  visual  representation  of  the  optimization 
perfonned  on  each  structure.  The  thesis  concludes  by  answering  the  questions  posed  in 
the  objective  section  along  with  recommendations  for  further  study.  The  appendix 
contains  samples  of  the  Matlab  code  used  to  perform  each  part  of  this  study. 

II.  Literature  Review 


2.1  Chapter  Overview 

There  is  a  scarcity  of  literature  available  on  vacuum  LTA  vehicles,  which  is  not 
surprising  considering  none  have  yet  been  built.  However,  each  of  the  components 
required  for  a  successful  vacuum  LTA  vehicle  have  been  extensively  studied  and 
published.  The  focus  of  this  section  is  to  introduce  LTA  vehicles,  their  structures, 
missions,  and  potential  areas  for  improvement. 

Although  a  vacuum  LTA  structure  has  apparently  never  been  built,  there  has  been 
extensive  research  into  the  buckling  of  pressure  vessels  under  external  pressure.  This 
section  examines  the  empirically  derived  fonnulas  for  failure  of  a  sphere  and  applies 
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them  to  the  problem  of  achieving  positive  buoyancy  in  air.  After  establishing  that  a  thin 
shelled  sphere  constructed  of  contemporary  materials  cannot  achieve  positive  buoyancy 
in  air  using  a  vacuum,  this  section  presents  two  possible  methods  of  geometrically 
stiffening  a  sphere  to  prevent  buckling  without  adding  too  much  weight.  These  methods 
are  an  isogrid  stiffened  sphere  and  a  geodesic  sphere. 

2.2  Historical  LTA  Vehicles 

LTA  vehicles  have  a  long  history  of  impressive  achievements  and  yet  have  been 
relegated  to  filling  niche  mission  profiles.  These  achievements  include  the  first  human 
flight,  accomplished  in  a  hot  air  balloon  by  Jacques-Etienne  Montgolfier  in  1782  [8,  p.  2] 
and  the  first  passenger  airline  Deutsche  Luftschiffahrts-AktiengeseUschaft  (DELAG)  in 
1909  [9],  In  1917  the  L-59  “Africa  Ship”  proved  the  feasibility  of  intercontinental 
zeppelin  travel  by  carrying  15  tons  of  cargo  and  22  persons  on  a  4,225  mile  flight  during 
a  supply  mission  to  German  East  Africa  [9],  The  qualities  of  LTA  vehicles  enabled  them 
to  accomplish  missions  unsuitable  for  other  vehicles,  but  also  historically  prevented  them 
from  competing  against  mainstream  methods  of  transportation.  Of  primary  importance  to 
this  study  is  airship  structure,  mission  profiles,  and  potential  for  improvement  on 
historical  designs. 

2.2.1  Airship  Structures 

Airship  structures  can  be  organized  into  three  primary  categories;  rigid,  semi¬ 
rigid,  and  non-rigid.  Each  of  these  structures  have  qualities  that  significantly  impact 
operational  capabilities. 
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2.2. 1.1  Rigid 

Rigid  structures  consist  of  a  rigid  frame  with  a  skin  stretched  over  the  outside  as 
shown  in  Figure  2  [10].  The  frame  distributes  loads  throughout  the  ship.  The  original 
Zeppelins  consisted  of  an  aluminum  frame  surrounded  by  a  penneable  cloth  skin. 
Hydrogen  lifting  gas  was  contained  within  animal  skin  ballonets  inside  the  structure. 
Weight/Buoyancy  (W/B)  for  these  structures  was  poor,  requiring  large  volume  to  achieve 
useful  designs.  Modem  material  and  construction  technologies  can  significantly  decrease 
the  weight  of  rigid  frames  and  skins  through  the  use  of  carbon  composites  and  advanced 
fabrics  and  membranes.  A  modem  example  of  a  rigid  airship  is  the  experimental 
Aeroscraft  created  by  the  Aeros  company  [11].  Rigid  structures  are  able  to  maintain  their 
shape  because  of  the  rigid  nature  of  their  structures.  This  enables  the  design  of  lifting 
bodies  as  well  as  the  ability  to  withstand  large  aerodynamic  loads  without  significant 
warping.  Rigid  structures  also  transmit  loads  throughout  the  structure  which  enables 
mounting  of  engines,  gondolas,  and  control  surfaces.  This  enables  improved 
controllability  by  side  mounted,  vectored  thmst  engines,  reduced  drag  by  internally 
placing  payloads,  and  optional  hybrid  lift  by  taking  advantage  of  a  lifting  body  shape. 
Hybrid  lift  is  the  use  of  aerodynamic  lift  in  addition  to  buoyant  lift.  Although 
aerodynamic  lift  incurs  significant  drag  in  the  case  of  a  LTA  vehicle,  it  allows  take  off 
and  landing  in  a  heavier  than  air  configuration  which  improves  ground  handling. 


9 


Figure  2:  USS  Shenandoah  Rigid  Construction  [9] 


Disadvantages  to  the  rigid  airship  design  include  increased  weight  for  the  frame 
[8,  p.  22]  and  poor  damage  tolerance  due  to  the  rigidity  of  the  structure  [9,  p.  150]. 
Whereas  a  non-rigid  structure  can  bend  without  damaging  the  structure  when  design 
bending  moments  are  exceeded,  a  rigid  structure  will  be  damaged  when  bending 
moments  are  exceeded. 

An  interesting  variation  of  the  rigid  airship  design  is  the  metal-clad  airship  as 
shown  in  Figure  3  [12].  The  US  Navy  Metal  Clad  Airship  (ZMC)-2  consisted  of  an 
aluminum  skin  with  aluminum  stringers  riveted  together.  The  ZMC-2  was  notable  for  its 
ability  to  retain  the  helium  lifting  gas  much  longer  than  contemporary  LTA  vehicles  [12]. 


This  was  due  to  the  low  penneability  of  the  aluminum  skin.  It  was  also  notable  for  its  low 
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fineness  ratio  (A),  which  is  defined  by  equation  (11)  and  is  the  ratio  of  length  to  diameter, 
which  was  chosen  to  add  strength  to  the  structure  [8,  p.  57].  The  low  A  strengthened  the 
ZMC-2  by  decreasing  bending  moments  and  increasing  the  resistance  to  bending  [13,  p. 
163],  [12].  There  are  published  equations  for  calculating  moments  on  an  airship  in  [13,  p. 
163],  but  a  simplified  explanation  as  to  why  low  A  decreases  moments  and  increases 
resistance  to  moments  can  be  made  by  picturing  the  LTA  vehicle  as  a  cylinder.  The 
longer  the  cylinder,  and  thus  the  higher  A,  the  higher  the  bending  moments  the  LTA 
vehicle  will  experience  simply  due  to  a  longer  moment  ann.  For  a  given  volume,  a  higher 
A  corresponds  to  a  lower  radius,  which  in  terms  of  a  cylinder  would  mean  less  resistance 
to  bending.  Also,  as  length  increases  and  volume  remains  the  same,  cross  sectional  area 
also  increases,  thus  increasing  loading  due  to  gusts.  A  low  A  decreases  stability  because 
of  the  lower  rotational  inertia  and  moment  arm  when  compared  to  a  higher  A,  requiring 
relatively  large  stabilizers  on  the  ZMC-2. 
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Figure  3:  ZMC-2  Metal  Clad  Airship  ca.  1930  |12] 


2.2. 1 .2  Semi-Rigid 

Semi-rigid  airships  consist  of  a  rigid  internal  frame  surrounded  by  an  inflated  skin 
that  maintains  shape  through  the  use  of  internal  pressure.  These  structures  carry  the 
advantages  of  both  the  rigid  and  non-rigid  designs.  The  rigid  internal  frame  distributes 
loads  throughout  the  ship  and  enables  freedom  to  mount  engines,  payloads,  and  control 
surfaces  anywhere  on  the  airship.  The  pressure  supported  skin  enables  a  smaller  internal 
frame  that  is  strength,  not  shape,  oriented,  thus  reducing  weight. 

A  modern  example  of  the  semi-rigid  airship  is  the  Zeppelin-NT,  which  has  an 

internal  frame  as  shown  in  Figure  4  [14].  The  Zeppelin-NT  rigid  frame  is  composed  of 

triangular  trusses  that  run  the  length  of  the  airship.  The  skin  touches  the  frame  at  specific 

points  but  maintains  its  shape  through  internal  pressure.  This  airship  can  conduct  ground 

handling  with  a  crew  of  three  people  using  vectored  thrust  [15].  Vectored  thrust  is  made 
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possible  by  the  side  mounting  of  its  engines,  which  is  made  possible  by  the  internal 
frame. 

2.2. 1.3  Non-Rigid 

Non-rigid  airships  consist  of  an  inflated  gas  envelope  that  maintains  shape 
through  internal  pressure.  These  structures  are  popularly  referred  to  as  “blimps”.  In  order 
to  withstand  mooring  and  aerodynamic  forces,  the  bows  of  non-rigid  airships  generally 
have  rigid  caps.  There  is  also  usually  a  rigid  structure  at  the  stem  for  mounting  the 
stabilizers.  The  gondola  generally  hangs  externally  from  a  canopy  that  distributes  forces 
over  a  large  portion  of  the  gas  envelope  as  shown  in  Figure  4  [14].  A  modem  example  of 
a  non-rigid  airship  is  the  1960s  era  Goodyear  Blimp  as  shown  in  Figure  4.  Non-rigid 
airships  do  not  require  significant  stmctures  and  therefore  can  achieve  low  W/B  and  high 
altitude  such  as  the  Lockheed  Martin  High  Altitude  Long  Endurance  Demonstrator.  They 
also  have  the  added  advantage  of  being  resilient.  If  the  design  bending  moment  of  a  non- 
rigid  airship  is  exceeded,  the  skin  may  temporarily  kink,  but  unlike  the  rigid  airship,  it 
will  not  be  damaged[13,  p.  164].  This  is  because  internal  gas  pressure  provides  stmctural 
support  as  opposed  to  a  rigid  stmcture. 
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Figure  4:  Semirigid,  Nonrigid,  Rigid  Airships  [14] 
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2.2.2  Missions 

The  original  mission  of  the  LTA  vehicle  was  aerial  reconnaissance.  In  1806 
Frenchman  Jean  Coutelle  was  lifted  450  meters  in  a  tethered  balloon  to  observe  enemy 
formations  and  movements  in  the  French  Revolution  [8,  p.  2].  Airships  were  used  in  the 
American  Civil  War,  the  Franco-Prussian  Wars,  World  War  I  and  World  War  II  [8,  p.  2], 
During  that  time  Airships  conducted  passenger  and  cargo  transportation,  reconnaissance, 
and  bombing.  As  heavier  than  air  flight  technology  matured,  LTA  vehicles  were  replaced 
by  heavier  than  air  vehicles  for  many  of  these  missions.  Today  LTA  vehicles,  to  include 
unmanned  airships  and  aerostats,  are  primarily  used  as  sensor  platforms  and 
communications  relays  [16].  Airships  in  particular  are  used  for  product  recognition, 
advertising,  tourism,  and  as  commercial  aerial  photography  platforms  [17].  The 
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endurance,  lift  capacity,  and  stability  of  LTA  vehicles  makes  them  well  suited  for  these 
missions,  but  there  are  aspects  of  LTA  vehicles  that  prevent  them  from  being  used  for 
cargo  or  passenger  transportation.  These  aspects  are  design  challenges  that  must  be 
overcome  if  the  LTA  vehicle  is  to  become  a  viable  competitor  in  the  transportation 
establishment. 

2.2.3  Design  Challenges 

The  first  design  challenge  of  LTA  vehicles  is  ground  handling,  which  historically 
requires  robust  hangars  and  mooring  capabilities.  A  related  design  challenge  is  buoyancy 
control  which  significantly  affects  mooring  and  taxiing.  Lifting  gases  pose  unique  supply 
and  operational  problems  that  must  be  addressed.  Efficiency,  as  it  relates  to  velocity  and 
fuel  requirements,  is  the  final  design  challenge. 

An  important  aspect  of  LTA  history  is  written  in  the  mishap  record.  The  mishap 
record  has  recurring  themes  that  shed  light  on  the  design  challenges  of  LTA  vehicles.  It  is 
the  duty  of  the  operator  to  recognize  these  themes  in  order  to  avoid  repeating  the 
mistakes  of  the  past,  and  the  duty  of  the  engineer  to  prevent  future  loss  of  life  and 
property  through  the  improvement  of  design.  Because  this  is  a  design  oriented  study,  it 
proved  fruitful  to  study  the  mishap  record.  The  mishap  study  gave  clear  answers  to 
several  of  the  questions  posed  in  the  objective  section  about  how  to  make  the  LTA 
vehicle  competitive  as  a  cargo  or  passenger  transportation  platform. 

A  review  of  mishap  records  for  LTA  vehicles  over  the  previous  century  revealed 

an  interesting  trend  as  shown  in  Figure  5  [9],  [18].  Mishap  causal  factors  for  LTA 

vehicles  are  shown  by  percent  of  total  reported  mishaps  for  both  hydrogen  and  helium 

airships  over  the  past  century.  It  is  interesting  to  note  that  pilot  error  accounted  for  far 
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less  mishaps  in  hydrogen  airships  than  in  helium  airships.  The  reason  for  this  trend  is 
possibly  twofold.  Pilots  of  hydrogen  airships  flew  at  a  time  when  airships  were  far  more 
prevalent,  and  training  was  much  more  robust  in  the  fonn  of  master  and  apprentice  [9], 
Pilots  of  helium  airships  were  flying  during  a  time  where  airships  were  one  of  a  kind, 
experimental  aircraft  with  inexperienced  crew,  resulting  in  a  higher  probability  of  pilot 
error  (with  the  exception  of  the  Goodyear  Blimps,  which  have  an  impressive  safety 
record  and  experience  pilots).  The  most  staggering  statistic  is  that  96%  of  reported 
hydrogen  airship  mishaps  involved  fire.  It  is  likely  that  there  were  many  unreported 
mishaps  and  that  those  mishaps  not  resulting  in  fire  were  often  not  important  enough  to 
be  included  in  the  historical  reports  used  for  this  study.  The  helium  airship  mishap 
statistics  come  from  modem  aircraft  mishap  reporting  requirements,  which  include  a  high 
number  of  minor  incidents.  This  skews  the  statistics  towards  making  reported  hydrogen 
airship  mishaps  seem  more  catastrophic  on  average  than  reported  helium  airship  mishaps. 
That  being  said,  the  number  of  hydrogen  airship  mishaps  resulting  in  death  is  a 
staggering  41%.  Even  more  devastating,  15%  of  Hydrogen  mishaps  resulted  in  100%  loss 
of  life  onboard.  Even  when  a  helium  airship  mishap  resulted  in  death,  it  was  usually  only 
one  or  two  people  with  the  vast  majority  of  those  involved  surviving  with  little  or  no 
injury.  In  addition  to  the  loss  of  life,  every  hydrogen  airship  mishap  involving  fire 
resulted  in  a  complete  destruction  of  the  airship,  which  represents  significant  monetary 
loss,  often  including  buildings  on  the  ground. 

Another  skew  to  the  statistics  is  that  many  of  the  mishaps  involving  maintenance 
did  not  involve  people  on  board  the  airships,  so  that  the  number  of  hydrogen  airship 


mishaps  not  resulting  in  death  or  injury  was  increased.  The  potential  for  destruction  by 
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hydrogen  explosions  and  fire  was  so  great  that  in  one  incident  four  airships  were 
destroyed  in  a  single  hangar  fire.  Goodyear  Blimps  stopped  using  hydrogen  as  a  lifting 
gas  after  the  Wingfoot  Air  Express  crashed  into  the  Chicago  Trust  and  Savings  Building 
killing  13  people  and  injuring  27  in  1919  [9].  The  US  government  stopped  using 
hydrogen  as  a  lifting  gas  after  the  US  Anny  airship  Roma  hit  power  lines,  ignited  and 
killed  34  out  of  45  people  on  board  in  1922  [9].  The  infamous  Hindenburg  was  actually 
designed  to  carry  helium,  but  due  to  the  US  policy  of  hoarding  helium  in  anticipation  of 
World  War  II,  the  US  refused  to  sell  helium  to  Nazi  Gennany  [9].  This  decision  resulted 
in  the  most  widely  publicized  airship  disaster.  Interestingly,  the  USS  Akron,  the  airship 
with  the  most  deadly  airship  crash  in  history,  was  a  helium  airship.  In  the  Akron  disaster, 
it  was  hypothermia  and  drowning  in  the  ocean  due  to  a  lack  of  life  boats  that  resulted  in 
the  deaths  of  73  out  of  76  crewmembers  [9], 
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Figure  5:  Causes  of  29  Helium  and  27  Hydrogen  Airship  Mishaps 
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The  study  of  mishaps  enables  the  airship  designer  to  avoid  past  mistakes  and 
mitigate  risk.  The  primary  cause  of  death  and  destruction  in  airship  history  is 
unequivocally  fire  as  a  result  of  using  hydrogen  as  a  lifting  gas.  The  simple  solution  to 
this  problem  is  to  use  an  inert  lifting  gas  such  as  helium.  Unfortunately,  if  helium  is 
unavailable  as  in  the  case  of  the  Hindenburg,  there  are  few  lifting  gas  alternatives.  This  is 
a  strong  argument  for  the  investigation  of  a  vacuum  in  place  of  a  lifting  gas. 

Weather,  pilot  error,  and  maintainer  error  are  all  mishap  causes  that  can  be 
mitigated  through  modern  aviation  practices.  This  leaves  ground  handling  and  buoyancy 
control  as  the  primary  safety  concerns  for  the  designer  of  LTA  vehicles.  These  two  issues 

are  actually  closely  related  and  can  both  be  mitigated  by  buoyancy  control. 

18 


2.2.4. 1  Ground  Handling  and  Buoyancy  Control 
The  problem  with  ground  handling  of  LTA  vehicles  is  that  they  are  large  and 
therefore  experience  great  forces  due  to  wind  gusts  that  can  push  them  into  obstacles, 
damage  their  structure  and  break  mooring  lines.  In  addition  to  large  wind  loads,  LTA 
vehicles  are  generally  close  to  neutral  buoyancy,  which  results  in  very  little  force 
anchoring  them  to  the  ground,  requiring  significant  mooring  measures.  Ground  handling 
of  the  large  airships  built  during  the  early  20th  century  required  crews  as  large  as  700 
men  [19,  p.  300]  to  hold  the  mooring  lines  as  shown  in  Figure  6,  [9].  These  operations 
were  so  labor  intensive  that  airships  had  to  land  or  take  off  in  the  early  morning  or  late 
evening  when  the  mooring  crews  were  off  work  [9].  The  problem  was  that  the  ground 
crews  were  not  full  time  airship  handlers,  they  had  full  time  jobs  and  were  part  time 
airship  handlers.  Many  airships  were  damaged  during  mooring  operations,  especially 
when  being  moved  in  and  out  of  hangars.  As  shown  in  Figure  5  about  a  quarter  of 
mishaps  occurred  during  ground  handling.  To  reduce  the  chance  of  mishaps,  anchor 
airships  for  servicing  between  flights,  and  decrease  manpower  requirements,  fixed  and 
portable  mooring  masts  were  constructed.  Portable  mooring  towers  were  originally  trains 
on  tracks  that  could  safely  move  the  airships  in  and  out  of  their  hangars.  Today’s  smaller 
LTA  vehicles  are  moored  by  trucks  with  towers  attached.  The  Zeppelin  NT,  through  the 
use  of  vectored  thrust  and  portable  mooring  towers,  can  moor  with  a  crew  as  small  as 
three  people  [15]. 
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Figure  6:  Ground  Handling  [9] 


In  order  for  LTA  vehicles  to  become  more  practical  they  would  require  the  ability 
to  land  at  small,  unimproved  airfields  without  mooring  facilities.  This  could  be 
accomplished  with  the  ability  to  land  and  taxi  heavier  than  air.  This  would  be  possible 
with  advanced  buoyancy  control  through  the  use  of  onboard  compressors  coupled  with 
hybrid  lift  for  landing  heavier  than  air.  If  the  buoyancy  of  a  LTA  vehicle  could  be 
precisely  and  rapidly  adjusted,  it  would  enable  vertical  take  offs  and  landings  without 
significant  ground  crew  involvement  or  mooring  equipment.  This  would  enable  use  of 
LTA  vehicles  at  unimproved  landing  sites,  greatly  increasing  their  usefulness.  This  point 
was  not  lost  on  Pasternak,  owner  of  the  Aeros  company,  who  is  currently  focusing  on 
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designing  the  next  generation  of  LTA  vehicles.  Aeros’  primary  technological 
improvement  over  traditional  designs  is  buoyancy  control  [20]. 

2.2. 4.2  Lifting  Gas 

All  methods  of  achieving  LTA  flight  to  date  have  relied  upon  lifting  gas. 
Common  lifting  gases  are  heated  air,  helium,  and  hydrogen.  Lifting  gases  serve  to 
displace  the  ambient  air  and  thus  achieve  positive  buoyancy.  LTA  flight  is  possible  due 
to  buoyancy,  which  is  defined  by  the  Archimedes  principal  which  states,  “the  buoyant 
force  on  a  submerged  object  is  equal  to  the  weight  of  the  fluid  displaced  by  the  object.” 
[8,  p.  14].  If  a  structure  will  float  in  air  the  weight  of  the  air  displaced  must  be  greater 
than  the  weight  of  the  structure  that  displaces  the  air.  The  weight  of  the  structure  will 
include  any  gas  contained  inside  the  structure,  so  that  a  structure  will  be  more  likely  to 
float  if  it  contains  a  gas  with  a  density  lower  than  the  air  it  displaces.  The  densities  of  the 
most  common  lifting  gases  are  plotted  in  Figure  7.  Hydrogen  has  the  lowest  density  of 
any  lifting  gas.  Helium  is  a  close  second,  while  heated  air  is  only  slightly  less  dense  than 
air  at  standard  atmospheric  temperatures.  In  the  case  of  an  absolute  vacuum  there  is  no 
lifting  gas  so  only  the  weight  of  the  structure  is  considered  to  detennine  whether  there  is 
positive  buoyancy. 
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Figure  7:  Comparison  Between  Lifting  Gases 


Elevation,  Thousands  of  Feet  MSL 


The  first  human  flight  was  conducted  in  a  LTA  vehicle  using  heated  air  as  a 
lifting  gas  [8,  p.  2],  Today  heated  air  is  only  used  as  a  lifting  gas  for  recreational  purposes 
due  to  its  relatively  low  buoyancy  and  high  energy  requirements.  The  high  energy 
requirement  stems  from  the  need  to  keep  adding  energy  to  the  lifting  gas  to  keep  the 
temperature  elevated. 

Hydrogen  has  the  lowest  density  of  any  lifting  gas  and  was  used  as  early  as  the 
civil  war  due  to  its  high  availability  [8,  p.  16].  Unfortunately  hydrogen  tends  to  produce 
violent  chemical  reactions  with  oxygen  in  the  presence  of  a  flame  or  spark,  which 


generally  results  in  catastrophic  mishaps  in  LTA  vehicles  as  shown  in  Figure  5,  [9],  For 
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this  reason  hydrogen  is  seldom  used  as  a  lifting  gas  today.  In  order  for  hydrogen  to  be  a 
viable  lifting  gas,  the  hydrogen  airship  would  have  to  be  designed  in  such  a  way  that  even 
in  a  mishap  the  hydrogen  would  not  be  ignited.  It  is  possible  that  an  airship  could  be 
designed  with  crash  resistant,  internal  ballonets  of  hydrogen  that  could  reduce  the 
chances  of  catastrophic  failure  to  an  acceptable  level.  Modem  buoyancy  control  would 
negate  the  need  to  valve  off  hydrogen  during  flight  which  would  prevent  dangerous 
mixtures  of  hydrogen  and  oxygen.  Non-flammable  membranes  could  also  be  used  to 
prevent  fire  from  spreading  from  an  engine  to  the  lifting  gas.  A  cost  benefit  analysis 
would  need  to  be  performed  to  determine  whether  these  additional  design  constraints 
would  make  hydrogen  a  reasonable  substitute  for  helium.  Regardless  of  how  many  safety 
measures  are  in  place,  hydrogen  always  has  the  potential  to  burn,  which  may  prohibit  it 
from  ever  being  a  viable  lifting  gas  where  human  life  is  at  risk. 

Helium  has  become  the  lifting  gas  of  choice  for  LTA  vehicles  because  it  is 

chemically  inert  and  has  the  second  lowest  density  of  any  gas  [8,  p.  16].  Problems  with 

helium  as  a  lifting  gas  are  primarily  the  small  size  of  the  helium  molecule,  which  tends  to 

pass  through  the  thin  membranes  used  to  enclose  it  [21,  p.  3],  as  well  as  mining  and 

shipping  costs.  These  problems  do  not  generally  pose  a  functional  obstacle  to  the  use  of 

helium  as  a  lifting  gas,  but  do  pose  economic  and  supply  challenges.  Furthennore,  the 

cost  of  helium  is  forecast  to  rise  significantly  in  the  next  few  years  due  to  increased 

global  demand  and  reduced  supply  [22].  Right  now  helium  is  widely  available  and  costs 

about  five  times  that  of  hydrogen  [23],  Helium  is  a  nonrenewable  resource  that  is  the 

product  of  radioactive  decay,  and  unlike  hydrogen,  can  not  be  manufactured.  At  current 

consumption  rates  it  appears  helium  will  continue  to  be  available  for  at  least  a  couple 
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more  decades,  however  costs  are  forecast  to  rise  significantly  in  the  immediate  future 
when  the  Bureau  of  Land  Management  stops  selling  helium  to  the  public  [23],  [22], 

2.2. 4. 3  Efficiency 

The  efficiency  of  any  vehicle  has  a  strong  influence  on  whether  or  not  it  is 
economically  feasible.  Even  if  the  handling  characteristics  of  LTA  vehicles  are  improved, 
the  velocity,  payload  and  efficiency  of  LTA  vehicles  will  have  to  be  comparable  to 
vehicles  with  similar  missions.  For  example,  if  the  LTA  vehicle  travels  the  same  velocity 
as  a  truck,  lands  at  the  same  distribution  centers  as  a  truck  ,  moves  the  same  payload  as  a 
truck,  but  uses  twice  as  much  fuel,  it  will  not  be  considered  a  viable  alternative.  Drag 
reduction,  weight  reduction,  and  increased  propulsion  efficiency  all  have  the  potential  to 
improve  LTA  vehicle  efficiency. 

Drag  reduction  can  be  accomplished  through  selection  of  optimal  fineness  ratio, 
Boundary  Layer  Control  (BLC),  and  a  clean  profile.  Fineness  ratio  (A)  is  defined  in 
equation  (1 1)  as  airship  Length/Maximum  Diameter  [8,  p.  44].  The  longer  a  ship,  the  less 
pressure  drag  it  experiences  for  a  given  volume,  but  the  more  skin  friction  drag  it  incurs. 
Therefore,  there  is  an  optimal  A  where  the  sum  of  skin  friction  drag  and  pressure  drag  is  a 
minimum.  A  A  =  4.62  was  found  to  result  in  the  lowest  total  drag  for  C-  Class  airships 
during  wind  tunnel  testing  [8,  p.  57],  but  the  ideal  A  will  change  with  various  design 
factors  such  as  aircraft  shape,  BFC  devices,  and  propulsion  location. 

BFC  can  be  used  to  reduce  pressure  drag  and  thereby  enable  a  lower  A  which  will 
in  turn  decrease  skin  friction  drag.  A  Griffith  airfoil  which  was  designed  to  achieve 
laminar  flow  through  BFC  is  shown  in  Figure  8  [24],  Laminar  flow  is  negligible  in 
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Reynolds  Numbers  (Re)  >  107  even  with  BLC,  but  BLC  in  the  case  of  an  airship  can 
prevent  flow  separation  and  thereby  decrease  pressure  drag.  Re  is  a  non-dimensional 
number  calculated  according  to  equation  (1),  where  Re  is  the  ratio  of  inertial  forces  to 
viscous  forces  [25,  p.  350].  An  airship  the  volume  of  the  Zeppelin  NT  at  8,297  m3  with  a 
A  —  2.5  would  be  47  meters  long.  If  it  was  travelling  at  a  um=36  m/s,  in  a  standard  sea 
level  atmosphere,  the  Re  would  be  about  1.15  ■  108,  meaning  a  Griffith  airfoil  would  not 
achieve  laminar  flow.  Goldschmied  [24]  conducted  wind  tunnel  testing  on  bodies  of 
revolution  based  on  the  Griffith  airfoil  where  he  showed  that  this  concept  for  a 
submerged  body  resulted  in  a  Cd=0.0162.  In  comparison,  the  Akron  Airship  model  had  a 
Cd  =0.0235  at  the  same  Reynolds  number.  In  addition  to  drag  reductions,  he  also  showed 
there  would  be  a  decrease  in  propulsive  power  required  by  45%  when  using  a  stem 
mounted  propeller.  His  design  consisted  of  a  combination  of  a  body  of  revolution 
designed  after  a  Griffith  airfoil  with  a  A  =  2.5,  which  maintained  a  positive  velocity 
gradient  along  the  body  of  the  airfoil  until  approximately  80%  chord,  followed  by  an 
abrupt  velocity  discontinuity.  A  slot  was  placed  just  aft  of  the  point  of  discontinuity 
which  prevented  the  flow  from  separating  using  suction  pumps.  These  pumps  then 
discharged  the  fluid  at  near  zero  velocity  into  a  stern  mounted  propeller,  thus  increasing 
the  propeller’s  efficiency  above  that  of  conventional  stern  mounted  propellers. 

n  Pair  '  ^oo  '  l  (1) 

Re  — - 

P 
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Where: 

l  —  airship  length 
itoo  =  freestream  velocity 
p  =  dynamic  viscosity 
Pair  =  air  density 


Figure  8:  Griffith  Airfoil  [24] 
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A  clean  profile  can  result  in  drag  coefficient  reductions  by  as  much  as  half.  The 
drag  coefficients  by  Goldschmied  [24]  were  for  a  clean  body  of  revolution  and  did  not 
include  a  gondola,  control  surfaces,  mooring  lines,  cables,  landing  gear,  antennas  etc... 
Through  the  use  of  modern  materials  most  of  those  objects  could  be  placed  inside  the 
envelope  of  the  airship,  significantly  decreasing  drag.  Control  surfaces,  however,  will 
most  likely  be  required  unless  control  can  be  achieved  through  a  creative  use  of  the 
suction  slots  and/or  vectored  thrust. 

Weight  reduction  is  important  in  any  aircraft  design.  In  the  case  of  the  airship  it  is 
useful  to  minimize  the  ratio  of  Weight/Buoyant  Force  (W/B).  A  W/B=l  would  be  a 
neutrally  buoyant  structure  that  neither  rises  nor  falls.  A  W/B<1  would  be  positively 
buoyant  and  float.  A  W/B>1  would  be  negatively  buoyant  and  sink.  The  lower  the  W/B 
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of  an  airship,  the  more  payload  it  can  carry  for  a  given  volume,  which  translates  in  to 
greater  transportation  efficiency.  Modem  materials  enable  the  construction  of  vehicles 
with  much  lower  W/B  than  in  the  past.  Modern  airships  often  achieve  W/B<0.3  [13,  p. 
418],  The  Zeppelin  NT  has  an  estimated  skin  plus  frame  W/B  S  0.18  [15]. 

Improved  propulsion  efficiency  can  be  accomplished  through  propeller  placement 
and  improved  engine  efficiency.  Stem  mounted  propellers  have  been  shown  to  reduce 
power  required  by  about  26%  by  placing  the  propeller  in  the  wake  of  the  airship  [24], 
There  is  an  additional  drag  savings  by  eliminating  the  need  for  fins  on  the  side  of  the 
airship  for  mounting  propellers.  A  stem  mounted  ducted  fan  or  propeller  in  conjunction 
with  a  BLC  suction  slot  on  a  Griffith  airfoil  was  predicted  to  reduce  power  required  by  as 
much  as  50%  [24].  The  tradeoff  with  stem  mounted  propulsion  is  a  loss  in 
maneuverability.  If  the  airship  gains  50%  efficiency,  but  loses  the  ability  to  land  at 
unimproved  airfields  it  may  not  be  a  viable  design.  Therefore  maneuverability  would 
need  to  be  considered  in  engine  placement.  A  possible  solution  is  to  use  propulsion 
similar  to  naval  ships,  which  often  have  both  bow  and  stem  vectored  thmst  [13,  p.  369]. 
By  also  placing  an  engine  near  the  bow  of  the  airship,  controllability  would  be  greater 
than  with  side  mounted  propulsion  due  to  increased  moment  arms. 

Engine  efficiency  can  possibly  be  improved  on  airships  through  the  use  of  high 
efficiency  engines  such  as  diesel  instead  of  high  specific  power  engines  such  as 
turboprops.  An  optimization  would  have  to  be  performed  which  took  into  account  the 
specific  mission  for  which  the  airship  was  being  designed  to  make  a  correct  judgment  on 
engine  type.  A  high  specific  power  engine  would  be  ideal  for  short  mission  durations,  but 
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a  high  efficiency  engine  might  be  preferable  for  long  missions  due  to  a  decreased  fuel 
weight. 

Efficiency  comparisons  between  vehicles  of  different  types  can  be  difficult.  A 
vehicle  designed  for  carrying  a  large  amount  of  cargo  can  be  expected  to  require  more 
power  than  a  vehicle  designed  for  carrying  a  few  passengers.  Also,  a  vehicle  designed  for 
high  velocities  can  be  expected  to  require  more  energy  than  one  designed  for  low 
velocity.  Gabrielli  and  von  Karman  defined  von  Kannan  efficiency  ( eVK )  according  to 
equation  (2)  [26].  Von  Karman  efficiency  enables  side  by  side  comparisons  of  different 
vehicles  by  calculating  efficiency  based  on  the  power  required  to  move  the  maximum 
payload  at  the  maximum  velocity  attainable  by  that  vehicle.  Note  that  eVK  is  a  non- 
dimensional  quantity  that  essentially  gives  a  ratio  of  power  required  to  useful  power 
output.  In  the  business  of  transporting  cargo  useful  work  would  be  defined  as  moving 
cargo,  and  useful  power  can  be  defined  as  moving  cargo  at  a  given  velocity.  Increasing 
velocity  is  desirable  for  some  missions,  but  there  is  an  associated  cost  of  lost  efficiency. 
This  relationship  between  efficiency  and  velocity  is  shown  in  Figure  9.  The  cost  of 
increasing  velocity  was  shown  by  the  Gabrielli-von  Kannan  limit  defined  by  equation 
(3).  A  plot  of  the  Gabrielli-von  Kannan  limit,  as  well  as  eVK  for  various  vehicles  is 
shown  in  Figure  9. 


eVK  — 


p 

rreq 

W  ■  um 


(2) 
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Von  Karman  Efficiency 


eVKLIMIT  ~  0.000175  ■  u 


(3) 


Where: 

evK  —  Von  Karman  Efficiency 
Preq  —  Power  Required 
u  =  velocity  mph 
um  =  Maximum  Velocity 
W  —  Useful  Lift 


Figure  9:  Von  Karman  Efficiency  of  BLC  Vehicles 


Figure  9  uses  data  for  various  vehicles  operating  with  their  typical  payload  and 


velocity  [26].  Larger  values  of  eVK  mean  less  efficient  modes  of  transportation  so  that  to 
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move  up  on  the  y  axis  is  to  be  less  efficient,  and  to  move  right  on  the  x  axis  is  to  achieve 
higher  velocity.  The  dotted  line  on  Figure  9  represents  the  Gabrielli-von  Kannan  Limit, 
(3).  This  limit  can  be  used  to  compare  the  relative  efficiency  of  different  vehicles.  A 
vehicle  that  falls  near  the  Gabrielli-von  Kannan  limit  can  be  considered  a  high  efficiency 
vehicle  in  comparison  to  other  vehicles  that  operate  at  the  same  velocity.  The  Gabrielli- 
von  Kannan  limit  is  based  on  viscous  drag  and  vehicle  data  available  in  1950.  This  limit 
can  be  broken  by  highly  efficient  modes  of  transportation  such  as  trains  and  tankers,  and 
so  should  be  thought  of  as  a  guideline  more  than  as  a  limit. 

The  dashed  curve  in  Figure  9  is  of  an  idealized  airship  with  an  internal  volume  of 
8,297  m3,  the  same  as  the  Zeppelin  NT  [15].  The  solid  line  is  the  same  idealized  airship 
design  with  an  internal  volume  of  77,956  m  .  The  larger  version  can  carry  a  payload  of 
13,617  kg  at  120  Knots  Indicated  Air  Speed  (KIAS)  at  10k  ft  for  1852  km  which  is  the 
range  for  the  C-27J  carrying  10,000  kg  of  cargo.  The  smaller  airship  represented  by  the 
dashed  curve  can  only  carry  1,361  kg  to  10k  ft  at  100  KIAS  for  1852  km.  This  is  a 
simplified  scenario  which  does  not  take  into  account  the  requirement  for  the  structure  to 
compensate  for  larger  engines  and  higher  aerodynamic  loads.  It  includes  increased  engine 
and  fuel  weight  according  to  equations  (6)  and  (7).  It  also  assumes  a  completely  clean 
aerodynamic  profile  where  the  landing  gear  and  control  car  are  shrouded  inside  the 
airship  to  reduce  drag.  The  curve  does  not  represent  the  performance  of  one  vehicle  at 
various  velocities  as  was  done  in  [26],  but  instead  the  performance  of  one  vehicle  at  its 
maximum  velocity  for  hundreds  of  different  sized  engines.  This  approach  was  taken  to 
explore  the  potential  design  space  of  a  theoretical  LTA  vehicle  in  order  to  determine  a 


suitable  maximum  velocity.  At  the  top  right  of  the  solid  and  dashed  curves  the  airship  has 
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an  engine  so  big  it  can  not  carry  any  payload  in  addition  to  the  required  fuel  and  engine, 
at  the  bottom  left  an  engine  so  small  it  can  only  propel  the  vehicle  at  10  kts.  The  solid  and 
dashed  curves  are  meant  to  represent  the  limit  of  what  can  be  expected  from  an  airship 
with  suction  type  BLC  that  is  completely  optimized  for  a  given  volume  for  weight,  drag, 
and  propulsion  efficiency.  This  shows  that  the  larger  the  airship,  the  more  efficient  it 
becomes.  It  also  shows  how  the  ideal  airship  compares  to  other  modes  of  transportation. 
Notice  that  the  von  Karman  efficiency  of  the  Zeppelin  NT  is  shown  in  Figure  9.  The 
potential  improvement  between  the  Zeppelin  NT  and  the  idealized  BLC  airship  with  a 
volume  equal  to  the  Zeppelin  NT  is  enormous.  The  actual  Zeppelin  NT  has  an  external 
gondola,  fin  mounted  propellers,  and  no  BLC.  It  cannot  be  reasonably  expected  that  an 
actual  vehicle  will  perform  as  well  as  the  solid  and  dashed  curves.  Instead  these  curves 
are  meant  to  represent  the  limit  of  conventional  airship  technology.  The  usefulness  of  the 
curves  is  that  they  show  the  advantage  of  designing  airships  for  low  velocities,  the  cost  of 
higher  velocities,  and  the  upper  limit  of  design  velocities.  This  limit  shows  that  above 
certain  design  velocities  an  increase  in  engine  power  only  serves  to  decrease  efficiency. 

Von  Karman  efficiency  was  calculated  for  the  idealized  airships  in  Figure  9  using 
Equations  (4)  -  (11).  The  assumptions  used  for  these  calculations  are  listed  in  Table  3. 
The  payload  weight  is  used  in  equation  (2)  and  calculated  using  equation  (4).  The 
payload  weight  is  the  buoyant  force  of  the  airship  minus  the  weight  of  the  structure, 
engine,  and  fuel.  Equations  (6)  [13,  pp.  225-226]  and  (7)  [25,  p.  233]  give  the  engine  and 
fuel  mass  which  are  based  on  the  power  requirements  of  equation  (8)  [25,  p.  238],  Engine 
mass  in  equation  (6)  is  based  on  average  specific  power  of  conventional  aircraft  turbine 


engines  after  they  have  been  installed  on  the  aircraft.  Power  is  a  function  of  volume  and 
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velocity.  The  plots  of  von  Karman  efficiency  vs.  velocity  were  made  possible  by  solving 
for  power  required  and  maximum  payload  weight  as  functions  of  velocity  for  a  given 
mission  range.  It  was  important  to  select  a  mission  range  because  this  determined  how 
much  fuel  needed  to  be  carried,  which  affected  the  payload  capacity.  Ranges  in  this  case 
were  based  on  the  C-27J  in  order  to  make  a  comparison  between  C-27J  capabilities  and 
conceptual  aircraft  designs.  The  drag  of  an  airship  as  shown  in  equation  (9)  is  based  on 
the  volumetric  drag  coefficient.  It  is  a  function  of  airship  volume  instead  of  frontal  area 
[27,  p.  4]  making  the  assumption  that  airships  have  the  classical  prolate  spheroid  shape. 
The  volume  of  an  airship,  which  makes  the  same  basic  shape  assumptions  as  equation 
(9),  is  shown  in  equation  (10)  and  is  a  function  of  fineness  ratio  A. 

W  =  Lb  —  g  ■  ( ms  +me  +  rrif )  (4) 
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Where: 

ap  —  propeller  area 

CD  —  Volumetric  Drag  Coefficient 

D  —  drag  force 

kee  —  1/engine  efficiency 

kpe  —  1/propeller  efficiency 

l  —  length 

Lb  —  buoyant  force 

me  =  engine  mass 

rrif  =  fuel  mass 

Px  =  power  required  to  maintain  translational  velocity 
r  =  radius 

V  —  volume  of  airship 

W  —  useful  load  assuming  W/B=0.3 

Z,  —  energy  density  of  fuel 

A  =  fineness  ratio 

Pair  =  density  of  air 

Pig  —  density  of  lifting  gas 


Table  3:  BLC  LTA  Vehicle  Assumptions 

ap  —  12.5  m2  Zeppelin  NT  sized  airship  [15] 

ap  =  56  m2  size  required  to  carry  10,000kg,  to  10k  ft,  1852  km 

CD  =  0.0169  based  on  hull  plus  tail  fins  [24] 

kee  —  0.33  33%  efficient  turboshaft  engines  [13,  p.  105] 
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kpe  =  ^-^7 assuming  stern  mounted  suction  slot  and  ducted  fan  [24] 

r  =  9.3,  m  Zeppelin  NT  sized  airship  [15] 
r  =  19.5,  m  size  required  to  carry  10,000kg,  to  10k  ft,  1852  km 
v  —  8,297  m3  volume  of  Zeppelin  NT  [15] 

v  —  77,956  m3  size  required  to  carry  10,000kg,  to  10k  ft,  1852  km 
ms  —  30%  of— assuming  W/B  =  0.3  [15] 

A  —  2.5  based  on  Griffith  Airfoil  [24] 

Pair  —  0.889  kg/m3  at  10k  ft  MSL  standard  atmosphere  [28] 

2.3  Thin  Shelled  Sphere 

LTA  structures  are  traditionally  composed  of  thin  membranes  filled  with  a  lifting 
gas.  If  the  density  of  the  lifting  gas,  ptg,  is  less  than  the  density  of  the  air  it  displaces, 
pair,  the  lifting  gas  is  said  to  have  positive  buoyancy.  If  the  buoyant  force,  Lb,  of  an 
object  as  calculated  by  equation  (5)  [8,  p.  15],  is  greater  than  the  weight  of  that  object, 
then  the  object  is  said  to  have  positive  buoyancy.  In  order  to  build  a  LTA  vehicle,  the 
weight  of  the  structure,  to  include  fuel  and  engines,  must  be  less  than  Lb  as  shown  in 
equation  (4),  otherwise  the  vehicle  will  not  have  positive  buoyancy.  The  dominant  weight 
in  equation  (4)  is  due  to  the  mass  of  the  structure,  ms.  The  design  of  a  vacuum  LTA 
structure  is  primarily  concerned  with  reducing  the  weight  of  the  structure  as  far  below  Lb 
as  possible.  In  order  to  determine  ms  the  geometric  and  material  properties  of  the 
structure  must  be  evaluated  such  as  geometric  stability  and  material  failure. 

A  vacuum  LTA  structure  is  one  in  which  a  partial  vacuum  is  used  instead  of  a 
lifting  gas.  Considering  equation  (5),  the  density  of  the  lifting  gas  would  be  determined 
by  the  ideal  gas  law,  (12)  [1]. 
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(12) 


Where: 

R  —  specific  gas  constant 

T=  temperature 

Pig  —  lifting  gas  pressure 

As  lifting  gas  pressure pig  approaches  zero,  lifting  gas  density  pLg  also  approaches 
zero.  The  pressure  of  the  lifting  gas  in  a  traditional  LTA  structure  is  generally  greater 
than  the  ambient  air  pressure  [13,  pp.  160-170].  This  results  in  skin  tension  which  adds 
stability  to  the  structure.  The  reverse  is  true  in  the  case  of  a  structure  that  contains  a 
vacuum.  In  this  case  the  external  pressure  on  the  structure  places  the  skin  under 
compression,  which  destabilizes  the  structure.  Structural  failure,  in  the  case  of  an 
externally  loaded  thin  shell,  can  occur  in  two  ways.  Buckling  is  the  first,  and  exceeding 
the  material  elastic  limit  is  the  second. 

Buckling  is  defined  as  an  inordinately  large  displacement  under  conditions  of 
unstable  equilibrium  [4,  p.  134].  A  thin  shell  in  compression  is  a  condition  of  unstable 
equilibrium  because  the  shell  has  the  potential  to  snap  through  locally  as  shown  in  Figure 
10  [29].  Figure  10  has  multiple  regions  of  local  buckling  instead  of  one  large  region  due 
to  artificial  stability  provided  by  a  solid  internal  mandrel.  After  buckling,  a  structure  will 
continue  to  displace  under  lower  loads  than  the  critical  load  at  which  buckling  first 
occurred.  If  buckling  occurs  there  will  be  two  primary  consequences.  One  is  material 
failure,  to  be  discussed  subsequently.  The  other  is  decreased  internal  volume,  which  will 
reduce  buoyancy.  This  buoyancy  reduction  can  be  sudden  if  buckling  occurs  and  could 


35 


result  in  instant  loss  of  positive  buoyancy.  Therefore  the  structure  must  avoid  buckling  if 
it  is  to  continue  to  accomplish  its  mission. 


Figure  10:  Buckled  Sphere  [29] 


The  critical  buckling  pressure  pcrit  of  a  homogeneous,  thin  walled,  sphere  under 
uniform  external  pressure,  is  shown  in  equation  (13)  [30,  p.  5]  which  is  based  on  classical 
buckling  theory.  In  practice  a  knockdown  factor  is  used  in  conjunction  with  equation  (13) 
to  match  experimental  results.  This  is  necessary  because  any  deviation  from  a  perfect 
sphere,  as  shown  in  Figure  1 1  [30,  p.  1 1],  results  in  a  lower  pcrit  than  theory  predicts. 
Figure  1 1  [30,  p.  11]  shows  the  knockdown  factor  for  the  buckling  pressure  of  a  prolate 
spheroid  plotted  against  fineness  ratio.  A  prolate  spheroid  is  a  body  of  revolution  similar 
in  shape  to  a  blimp  in  which  one  axis  of  the  spheroid  is  longer  than  the  other,  and  the 
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long  axis  is  the  axis  of  revolution.  An  absolutely  perfect  sphere  would  match  classical 
buckling  theory,  but  even  minor  deviations  cause  a  rapid  drop  off  in  pcrit  as  shown  in 
Figure  11.  Classical  buckling  theory  was  also  used  in  equation  (14)  [4,  p.  1292],  but  with 
the  addition  of  an  empirical  knockdown  factor  of  about  0.7  that  gives  a  more 
conservative  prediction  of  pcrit .  A  knockdown  factor  of  0.7  was  recommended  by 
Krenzke  [3 1,  p.  14]  after  experimental  analysis  of  near  perfect  spherical  caps  gave  a 
knockdown  factor  of  0.73  and  0.9  of  that  predicted  by  classical  theory  [3 1,  p.  9]. 


Pcrit  2  ' 


E 

V3(l-v2) 


(13) 


Pcrit  0.8  ■ 


(14) 


Where: 

E  =  Young’s  Modulus  of  Elasticity 
tskins  —  sphere  skin  thickness 
Pcrit  —  critical  pressure 
rs  =  sphere  radius 
v  =  Poisson’s  Ratio 
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Figure  11:  Prolate  Sphere  |30] 


B 

The  second  way  in  which  structural  failure  can  occur  is  when  the  elastic  limit  of  a 
material  is  exceeded.  In  order  to  ensure  a  sphere  will  not  collapse  due  to  material  failure, 
equation  (15)  [4,  p.  124]  must  also  be  satisfied.  This  was  derived  by  setting  the 
compressive  yield  strength  of  the  sphere  equal  to  the  load  divided  by  the  cross  sectional 
area  as  shown  in  equation  (16). 
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Where: 

ocirc  —  cross  sectional  area  of  great  circle  of  sphere 
cl  skin  —  cross  sectional  area  of  sphere  skin 
oCy  =  material  compressive  yield  strength 

A  sphere  that  can  achieve  positive  buoyancy  using  a  vacuum  must  have  a  mass 
less  than  the  air  it  displaces,  and  satisfy  the  thickness  and  material  property  requirements 
of  equations  (14)  and  (15).  The  mass  of  a  spherical,  homogeneous  structure  can  be 
detennined  by  multiplying  the  density  of  the  skin  by  the  volume  of  its  skin.  Equation  (19) 
gives  the  mass  of  a  spherical  structure  and  is  the  area  of  a  sphere  multiplied  by  its 
thickness  and  density. 


Uls  Ps  '  tskin  '  ‘A 
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(20) 


as  —  4  ■  n  ■  r2 


Where: 

as  —  surface  area  of  spehre 
ms=  mass  of  structure 
tskin  =  skin  thickness 
ps  =  density  of  sphere  wall 

The  volume  of  a  sphere  is  shown  in  equation  (21)  [5],  Inequality  (22)  is  based  on  a 
combination  of  equations  (14),  (19),  and  (21),  by  setting  the  W/B<1.  This  gives  the 
required  material  properties  for  a  homogeneous,  thin  walled  sphere  to  achieve  positive 
buoyancy  in  air. 
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By  assuming  standard  atmosphere  at  sea  level  and  v  =  0.3,  inequality  (22)  reduces  to 
inequality  (23)  [32]  which  can  be  used  to  detennine  whether  a  material  would  be  capable 
of  achieving  positive  buoyancy  using  a  vacuum  at  sea  level.  The  value  in  inequality  (23) 
is  stiffness  over  density  squared,  which  shows  that  the  desired  material  properties  for  a 
thin  shelled  sphere  are  a  low  density  and  a  high  stiffness.  The  same  approach  was  taken 
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in  [5]  where  the  value  of  inequality  (23)  was  found  to  be  4.5  ■  105.  This  is  less 
conservative  than  the  value  in  inequality  (23)  due  to  the  use  of  the  classical  buckling 
theory  without  a  knockdown  factor  as  shown  in  equation  (13). 

A  survey  of  published  material  properties  indicated  that  there  was  no  material 

available  that  meets  the  requirements  of  inequality  (23).  The  material  properties  of  two 

materials,  where  the  fourth  column  is  the  left  side  of  inequality  (23),  are  shown  in  Table 

4.  The  pitch  based  fiber  represents  the  highest  specific  stiffness  composite  material,  and 

beryllium  the  highest  specific  stiffness  isotropic  material  available  on  open  source 

material  data  bases  such  as  Matweb  [33].  As  Table  4  shows,  neither  material  achieves  a 

specific  modulus  high  enough  to  satisfy  the  requirements  of  inequality  (23).  There  are 

two  rows  devoted  to  Dupont  E-130-X  Pitch  Based  UHM  Carbon  Fiber.  The  first  shows 

the  longitudinal  stiffness  for  unidirectional  fiber.  Isotropic  properties  are  assumed  for 

equation  (14),  so  the  fiber  stiffness  cannot  be  used  directly.  A  [(0/±45/90)]s  ply  has 

quasi-isotropic  properties  but  the  effective  stiffness  is  about  38%  of  the  unidirectional 

composite  axial  stiffness  [34,  p.  143].  This  value  was  calculated  using  laminate  theory 

which  was  beyond  the  scope  of  this  study.  However,  a  simple  estimation  approach  is 

often  used  for  aerospace  applications  that  gives  nearly  identical  results  to  laminate 

theory.  This  method  uses  the  rule  of  mixtures  and  reinforcing  efficiency  as  described  in 

[35].  The  rule  of  mixtures  states  that  the  stiffness  of  the  composite  material  is  equal  to  the 

combined  stiffness  of  the  fibers  and  matrix  times  their  respective  volume  fractions  as 

shown  in  equation  (24).  Volume  fraction  is  the  fraction  of  the  total  volume  of  the 

composite  composed  of  either  fiber  (Vf)  or  matrix  (Vm).  In  the  case  of  carbon  epoxy  the 

carbon  fibers  can  make  up  as  much  as  60%  of  the  composite  volume,  resulting  in  a 
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Vf  —  0.6  and  a  Vm  —  0.4.  The  fiber  reinforcing  efficiency,  rjg,  is  based  on  the  angle  of 
fibers  and  their  volumetric  portion  of  total  fiber  content  as  shown  in  equation  (25).  For  a 
[(0/±45/90)]s  ply,  which  gives  rig  =  0.375  according  to  equation  (25),  with  fiber 
stiffness  Ef  =  895  GPa,  and  matrix  stiffness  Em  =  3.5  GPa  at  Vf  =  0.6,  equation  (29) 
gives  a  composite  stiffness  of  Ec  —  0.375  ■  895  ■  0.6  +  3.5  ■  0.4  =  202  GPa.  This  value 
is  reflected  in  the  second  row  of  Table  4.  As  a  point  of  comparison  to  laminate  theory  as 
used  in  [34],  which  gave  a  quasi-isotropic  stiffness  of  38.2%  of  the  unidirectional 
stiffness,  the  quasi-isotrpoic  stiffness  calculated  using  equations  (24)  and  (25)  of  202  G 
Pa  was  27.5%  of  the  uniderectional  composite  stiffness  of  538  GPa. 

Ec  =  VeEfVf  +  EmVm  (24) 


„=Y  CLn  ■  COS4  9  (25) 


Where: 

an  —  proportion  of  total  fiber  content 

Ec=  composite  stiffness 

Ef=  fiber  stiffness  in  fiber  direction 

Em=  matrix  stiffness 

Vm  —  volume  fraction  of  matrix 

Vf  —  volume  fraction  of  fiber 

r]g  —  Krenchel  efficiency  factor 

9  —  angle  of  fibers 
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Table  4:  Material  Properties-Monocoque  Shell 


Material  from 
Matweb 

E,  Pa 

Ps? 

kg_ 

m3 

E 

pi 

Smooth 

Sphere 

W/B 

Comments 

Dupont  E-130-X 
Pitch  Based  UHM 
Carbon  Fiber 

895 ■ 109 

2130 

[(0/±45/90)2]s 

202 ■ 109 

1522 

8.7-  104 

3.5 

Composite 

Beryllium  SR-200 

303 ■ 109 

1840 

8.9-  104 

3.6 

Isotropic 

metal 

The  [(0/+45/90)]s  ply  was  chosen  because  it  is  quasi-isotropic.  This  means  that  in 
the  plane  of  the  laminate  it  gives  equal  material  properties  in  every  direction.  The 
stiffness  of  composites  made  of  Dupont  E-130-X  carbon  fiber  for  various  angles  with 
respect  to  the  axial  direction  as  calculated  by  equations  (29)  and  (30)  are  shown  in  Figure 
12.  These  curves  show  that  a  unidirectional  laminate  is  very  stiff  in  the  axial  direction, 
but  the  quickness  quickly  drops  off  at  other  orientations.  If  a  unidirectional  laminate,  or 
even  a  [0/90] s  ply  were  used  to  construct  a  sphere  there  would  be  directions  of  high 
stiffness  and  directions  of  low  stiffness  that  would  invalidate  the  assumptions  of  equation 
(14). 


43 


Figure  12:  Laminate  Stiffness  as  a  Function  of  Orientation 
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Spheres  under  external  pressure  may  fail  at  much  less  than  the  critical  pressure  of 
a  perfect  sphere,  or  even  less  than  the  knockdown  factor  of  0.7  [3 1,  p.  1].  Geometric 
stiffeners  and  increased  thickness  help  reduce  this  imperfection  sensitivity  [29], 
Considering  that  no  materials  can  enable  a  thin  shelled  sphere  to  achieve  positive 
buoyancy,  and  that  even  what  can  be  achieved  is  highly  susceptible  to  small 
imperfections,  another  method  is  required  to  accomplish  the  goals  of  this  thesis.  For  this 
reason  geometry  became  the  dominant  factor  in  the  search  for  a  structure  that  could 
achieve  positive  buoyancy  in  air  using  a  vacuum. 
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2.4  Isogrid 


Grid  stiffened  structures  consist  of  a  thin  skin,  as  in  the  thin  shelled  sphere 
evaluated  above,  with  a  grid  of  stiffeners  attached  to  the  skin.  If  the  grid  is  arranged  in 
such  a  way  that  the  stiffening  ribs  transmit  forces  equally  in  all  directions  tangent  to  the 
skin,  it  is  called  an  isogrid.  An  isogrid  structure  can  be  treated  as  an  isotropic  material 
which  enables  the  use  of  isotropic  equations  such  as  equation  (14)  [36,  p.  2.0.001].  The 
grid  is  generally  arranged  in  equilateral  triangles.  When  applied  to  a  spherical  surface 
they  can  be  arranged  in  the  fonn  of  a  geodesic  sphere  as  shown  in  Figure  13. 

Figure  13:  Grid  Stiffened  Sphere, /=10 


2.5  Geodesic  Sphere 

Geodesic  spheres  are  geometric  shapes  that  approximate  a  sphere  using  straight 

lines  along  the  great  circles  of  a  sphere  as  shown  in  Figure  13.  Buckminster  Fuller  was 

granted  a  patent  for  his  geodesic  dome  design  on  12  Dec  195 1  [37].  In  his  patent,  Fuller 
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described  a  geodesic  sphere  based  on  a  20  sided  polyhedron  called  an  icosahedron.  The 
icosahedron,  shown  in  Figure  14,  is  a  polyhedron  composed  of  20  equilateral  triangles  in 
which  each  vertex  lies  on  the  surface  of  an  imaginary  sphere  (shown  as  dashed  circle). 
Each  face  of  the  icosahedron  is  referred  to  as  a  major  triangle.  To  form  more  complex 
geometric  structures  that  more  closely  approximate  a  sphere,  the  major  triangles  are 
subdivided  into  smaller  triangles.  The  vertices  of  the  minor  triangles  that  subdivide  the 
major  triangles  also  lie  on  the  surface  of  the  sphere  that  circumscribes  the  vertices  of  the 
original  icosahedron.  The  number  of  divisions  along  one  edge  of  each  major  triangle  is 
referred  to  as  geometric  frequency  if).  A  geodesic  sphere  with  /=2  is  one  in  which  the 
edges  of  each  major  triangle  are  divided  into  two  as  shown  in  Figure  15.  An  increase  in/ 
results  in  increasingly  more  complex  structures.  Note  that  increasing/beyond  1  results  in 
a  geodesic  structure  composed  entirely  of  irregular  triangles,  whereas  the  original 
icosahedron  is  entirely  composed  of  equal  sized  equilateral  triangles.  This  point  becomes 
important  when  evaluating  loading  of  the  geodesic  sphere,  because  there  is  perfect 
symmetry  in  an  icosahedron,  but  not  on  more  complex  geodesic  spheres. 
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Figure  14:  Geodesic  Sphere /=  1  (Icosahedron) 


Figure  15:  Geodesic  Sphere /=  2 


A  Matlab  program  titled  LTAgeod.m,  Appendix  A.4,  was  used  to  create  Figure 


13,  Figure  14,  and  Figure  15.  This  program  generated  the  geodesic  sphere  geometry  by 
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starting  with  an  icosahedron  and  subdividing  each  major  triangle  according  to  f.  All 
coordinates  were  generated  in  spherical  coordinates,  but  later  converted  to  cartesian  for 
use  in  a  Finite  Element  Analysis  (FEA)  program.  Spherical  coordinates  are  defined  in  the 
Matlab  sph2cart.m  help  file  according  to  Figure  16  [38],  and  are  composed  of  6,  <p,  and  r. 
A  point  in  spherical  coordinates  is  defined  by  a  vector  that  has  length  r,  and  direction 
defined  by  6  and  <p.  The  angle  6  is  measured  counterclockwise  in  the  xy  cartesian  plane 
from  the  x  axis  to  the  position  vector.  The  angle  <p  is  measured  between  the  position 
vector  and  the  xy  cartesian  plane  as  shown  in  Figure  16. 


Figure  16:  Spherical  to  Cartesian  Conversion  in  Matlab  |38] 


x  =  r  .*  cos(phi)  .*  cos(theta) 
y  =  r  .*  cos(phi)  .*  sin(theta) 
z  =  r  .*  sin(phi) 


The  spherical  coordinates  for  the  vertices  of  an  icosahedron  with  center  at  the 
origin  are  listed  in  Table  5.  Every  vertex  of  an  icosahdron  lies  on  the  surface  of  an 
imaginary  sphere.  By  placing  the  center  of  the  sphere  at  the  origin,  the  coordinates  of 
every  point  on  the  surface  of  the  sphere,  to  include  the  vertices  of  the  icosahedron,  are  at 
r  =  rs. 
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Table  5:  Icosahedron  Vertex  Locations 

Point 

6  (radians) 

<p  (radians) 

r  (meters) 

1 

0 

n/2 

2 

0 

0.464 

r. 

3 

2?r/5 

0.464 

4 

47t/5 

0.464 

5 

6n/S 

0.464 

r 

6 

8n/S 

0.464 

fir 

7 

tt/5 

-0.464 

r. 

8 

3tt/5 

-0.464 

r. 

9 

n 

-0.464 

10 

7n/S 

-0.464 

fir 

11 

9n/5 

-0.464 

r. 

12 

0 

-  7t/2 

rs 

The  coordinates  in  Table  5  were  derived  with  the  knowledge  that  an  icosahedron 
is  composed  of  20  equilateral  triangles  with  vertices  circumscribed  by  a  sphere.  The 
center  of  the  icosahedron  was  placed  at  the  origin  with  the  z  axis  going  through  the  top 
and  bottom  vertices.  Thus  the  top  and  bottom  vertices  are  defined  by  0=±  -  and  6  —  0. 
The  rest  of  the  vertices  lie  on  one  of  two  planes  equidistant  from  the  xy  cartesian  plane, 
each  plane  containing  5  points  of  constant  <p  that  make  up  a  pentagon.  The  top  plane  is 
positive  0,  the  bottom  plane  is  negative  <p.  Therefore  <p  only  has  to  be  calculated  once. 
The  complete  mathematical  derivation  for  an  icosahedron  was  beyond  the  scope  of  this 
thesis,  however  starting  from  the  relationship  between  the  length  of  the  edge  of  a  major 
triangle  (a),  and  the  radius  of  the  circle  (rs)  as  shown  in  (26)  [39],  the  derivation  of  <p  was 
fairly  straight  forward.  By  cutting  an  icosahedron  in  half  using  a  plane  on  the  z  axis  that 
passed  through  a  vertex  on  the  top  pentagon  as  shown  in  Figure  17,  the  angle  <p  for  the 
points  on  the  top  and  bottom  pentagon  were  determined.  The  two  dimensional  view  of 
the  cut  icosahedron  is  shown  in  Figure  18.  The  angle  <p  can  be  found  using  the  triangle 
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formed  between  two  vertices  of  the  icosahedron  and  the  origin.  This  triangle  has  two  legs 
of  length  rs  and  one  leg  of  length  a  as  shown  in  Figure  18.  The  angle  opposite  the  leg  of 
length  a  was  detennined  by  the  law  of  cosines  as  shown  in  equation  (27)  [40],  and  came 

out  to  be  cos-1  s  1.107.  Therefore,  <p  =  ^  —  1.107  =  0.464. 


r  = 


27T 


a  ■  sin- 


(26) 


legf  —  leg |  +  legl  —  2  ■  leg2  ■  leg3  ■  cos (angle^ 


(27) 
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Figure  17:  Icosahedron-Cutting 


Figure  18:  Icosahedron-Cut 


z  axis 

A 


xy  plane 
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The  angle  0  was  easier  to  calculate.  Looking  at  the  icosahedron  from  top  down,  as 
in  Figure  24,  it  is  apparent  that  each  pentagon  divides  the  circle  into  five  equal  angles. 
Therefore,  6  comes  in  multiples  of  2  -n/S.  The  first  vertex  on  the  top  pentagon  is  placed 
on  the  x  axis  so  that  0  =  0.  The  following  vertices  on  the  top  pentagon  are  at  6  —  2  ■  n/ 
5,0  =  4  ■  7t/5,  0  —  6  ■  n/S  and  so  on  for  each  of  the  five  vertices  on  the  top  pentagon. 
The  top  and  bottom  pentagon  are  offset  from  each  other  so  that  each  angle  on  the  bottom 
pentagon  bisects  the  two  closest  angles  on  the  top  pentagon  as  shown  in  Figure  19. 
Therefore  0  on  the  bottom  pentagon  starts  at  0  —  n/S  and  increases  by  increments  of 
2  ■  7r/5  for  each  subsequent  vertex. 


Figure  19:  Icosahedron-Top  View 


a  y  axis 


_  _>  x  axis 


The  number  of  vertices  ( nv ),  edges  (ne),  and  triangular  faces  ( nt )  of  a  geodesic 
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structure  as  a  function  of f  are  calculated  using  equations  (28)  -  (30)  [41].  An 
icosahedron  has  an^=l,  which  gives  12  vertices,  30  edges,  and  20  triangles  using 
equations  (28)  -  (30).  These  equations  show  that  as/increases,  the  complexity  of  a  sphere 
increases  significantly.  The  number  of  vertices,  edges,  and  triangles  were  used  to  check 
the  program  LTAgeod.m,  Appendix  A. 4  for  accuracy,  and  for  weight  calculations  when 
solving  for  W/B. 


nv  =  10  ■  f2  +  2 


(28) 


nt  =  20  ■  f2 


(29) 


ne  =  30  •  / 2 


(30) 


2.6  Summary 

A  literature  review  was  conducted  that  covered  each  of  the  elements  required  to 
evaluate  a  vacuum  LTA  vehicle.  These  elements  included  general  LTA  vehicle  aspects  as 
they  pertained  to  the  problem  of  constructing  a  vacuum  LTA  vehicle.  These  included 
historical  structures,  missions,  and  potential  areas  for  improving  on  historical  designs. 
Thin  shelled  spheres  were  examined  which  showed  the  critical  material  properties  for 
building  a  vacuum  LTA  structure,  as  well  as  demonstrated  the  need  for  geometric 


53 


methods  of  preventing  buckling.  Isogrid  stiffened  spheres  and  geodesic  geometry  were 
introduced  as  possible  methods  of  preventing  buckling  without  adding  significant  weight. 


III.  Methodology 

3.1  Chapter  Overview 

This  chapter  introduces  the  rotating  cylinder  LTA  vehicle  concept  and  the 
methods  used  to  evaluate  its  feasibility.  These  include  buoyancy,  weight,  aerodynamics, 
power  requirements,  and  engine  and  fuel  weight.  Also  explained  are  the  equations  from 
the  NASA  Isogrid  Handbook  [36]  and  the  equations  and  methods  behind  the  programs 
used  to  optimize  the  geodesic  sphere. 

3.2  Rotating  Cylinder  Vacuum  LTA  Vehicle  Concept 

A  preliminary  investigation  into  the  feasibility  of  constructing  a  thin  shelled, 
unreinforced  structure  that  achieved  positive  buoyancy  in  air  using  a  vacuum,  showed 
that  there  was  no  material  in  existence  by  which  such  a  structure  could  be  fonned. 
Therefore,  it  became  apparent  that  the  structure  would  need  to  be  stiffened  geometrically 
or  by  some  other  means.  One  method  of  stiffening  the  structure  was  inspired  by 
helicopter  rotor  blades.  Helicopter  blades  are  highly  flexible  and  free  to  flap  vertically 
due  to  a  pinned-free  boundary  condition,  but  due  to  high  rotational  velocity  they  become 
resistant  to  vertical  movement  and  can  support  great  loads.  This  same  concept  was 
applied  to  the  vacuum  structure  by  investigating  long  cylinders  rotated  about  their  axis  of 
symmetry.  If  these  cylinders  were  rotated  at  a  sufficient  velocity,  the  centripetal  force 
exerted  by  the  skin  would  produce  a  hoop  tensile  stress  that  would  counteract  the  hoop 
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compressive  stress  in  the  skin  produced  by  atmospheric  pressure.  To  put  this  simply, 
resisting  atmospheric  pressure  pushes  the  skin  towards  the  center  while  inertia  pulls  the 
skin  away  from  the  center.  The  assumption  was  made  that  dynamic  effects  would  not 
interfere  with  the  balance  between  inertial  and  pressure  forces.  This  assumption  would 
need  to  be  verified  through  simulation  or  experiment,  but  it  allowed  energetic  evaluation 
of  the  rotating  cylinder  idea.  If  the  rotating  cylinder  proved  energetically  feasible,  the 
next  logical  step  would  be  to  evaluate  the  dynamic  effects  of  rotating  a  long,  slender 
cylinder.  Assuming  perfect  balance  between  all  radial  forces,  energy  requirements  were 
evaluated  based  on  drag  torque,  bearing  friction,  turbine  engine  efficiency,  propeller 
efficiency,  and  power  required  to  overcome  translational  drag. 

A  theoretical  vehicle,  Figure  20,  was  proposed  that  consisted  of  two  counter¬ 
rotating  cylinders  mounted  vertically  above  a  gondola.  The  gondola  would  have  two 
propellers  capable  of  thrust  vectoring.  At  the  top  of  the  counter-rotating  cylinders  would 
be  a  third  propeller  which  would  be  used  to  trim  the  vehicle  to  keep  the  rotating  cylinders 
perpendicular  to  the  airflow.  The  buoyancy,  weight,  aerodynamic,  and  power 
characteristics  of  this  vehicle  were  evaluated  to  determine  whether  it  was  a  feasible 
concept.  This  vehicle  concept  was  called  the  Rotating  Cylinder,  Vacuum  LTA  Vehicle 
(RCVLTAV).  Each  cylinder  was  proposed  to  be  a  smooth,  thin  shelled  structure  with 
spherical  end-caps  similar  to  commercial  pressure  vessels  designed  to  contain 
compressed  gas. 
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Figure  20:  Rotating  Cylinder  Vehicle 


3.2.1  Buoyancy  and  Weight 

The  buoyant  force  of  the  theoretical  structure  was  determined  using  equation  (3 1) 
by  calculating  the  volume  of  each  cylinder  with  equation  (32)  [42],  and  the  volume  of  the 
spherical  end-caps  with  equation  (33)  [42].  Both  cylinders  were  designed  to  be  of  equal 
dimensions.  The  buoyant  force  of  a  cylinder  with  spherical  end-caps  is  the  mass  of  air 
displaced  by  the  cylinder  and  end-caps  as  shown  in  equation  (31). 

Lbc  =  9  ‘  Pair  ■  Oc  +  Vs)  (31) 
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(32) 


vc  =  n-lc-  rc2 


v,  =  —  -n-r} 

s  3  c 

Lbc  —  buoyant  force  of  cylinder  with  spherical  endcaps 
uc=eylinder  volume 
us=spherical  endcap  volume 


(33) 


After  calculating  the  buoyancy  of  a  given  set  of  cylinders,  the  thickness  of  the 
cylinder  walls  was  determined  based  on  equation  (34)  which  was  derived  from  equations 
(35)  -  (37)  [4].  These  equations  prevent  buckling  of  the  cylinder  due  to  the  axial  load  on 
the  cylinder  caused  by  atmospheric  pressure  on  the  spherical  end-caps  as  depicted  in 
Figure  2 1 .  The  axial  force  on  a  cylinder  due  to  pressure  on  the  spherical  end-caps  minus 
the  buoyant  force  is  shown  in  equation  (35).  The  reason  the  buoyant  force  is  subtracted 
from  the  axial  load  is  that  the  payload  of  the  cylinder  hangs  from  one  end,  thus  placing 
the  cylinder  in  tension  by  an  amount  equal  to  the  buoyant  force.  Equation  (36)  is  the 
equation  for  the  buckling  end  load  of  a  hollow,  cylindrical  column  with  free-free 
boundary  conditions.  Using  the  end-load  calculated  in  equation  (35),  the  required 
thickness  for  the  cylinder  was  found  using  equations  (36)  and  (37)  [4,  p.  1296].  The 
symbolic  solution  to  these  equations  is  shown  as  equation  (34).  The  compressive  yield 
stress  of  the  cylinder  material  was  compared  to  equation  (38)  [4],  Cylinder  skin 
thickness,  tskirLc,  was  chosen  to  satisfy  both  equations  (34)  and  (38).  End-cap  thickness, 
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t Skms,  was  chosen  to  satisfy  equation  (14).  Material  properties  were  based  on  UHM 
carbon  epoxy  as  shown  in  Table  7. 


Lskinr  ~ 


1 64  ■  rc6  4  ■  Fa2  ■  Zc4  2  ■  Fa  ■  lc 

'  r  -w-.  n  n 


27  7T6  ■  E2  ■  rc2  n3  ■  E  ■  rc 


(34) 


64 -rc6  4  -Faz-lc  2  •  Fa  •  lc 
+  I  -  \  + 


27  n6  ■  E2  ■  rc2  n3-E-rc 


Fa  ~  Patm  m  ft  m  rc  L i 


(35) 


n2  ■  E  ■  I 


(36) 


I  —  n 


Lskin 


T  t skinc  '  Vc 


(37) 


On  = 


a  2  ■  n  ■  rr  ■ 


skinr 


(38) 


Where: 

Ea  —  axial  force  on  cylinder 

rc  —  cylinder  radius 

vc  =  cylinder  volume 

mc  —  cylinder  mass 

lc  —  cylinder  length 

tskinc  —  cylinder  skin  thickness 

tskins  —  sphere  skin  thickness 

oa  —  axial  stress 
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Figure  21:  Rotating  Cylinder  Loads 


After  determining  tskiUc  and  tskills ,  it  was  possible  to  calculate  the  mass  of  each 
cylinder  using  equation  (39)  where  the  area  of  the  cylinder  and  spherical  end-caps  were 
calculated  with  equations  (40)  and  (41)  [42].  The  assumption  was  made  in  equation  (39) 
that  there  was  an  absolute  vacuum  inside  of  each  cylinder.  This  assumption  provided 
reasonable  results,  even  though  an  absolute  vacuum  can  not  be  practically  achieved.  A 
very  low  pressure,  such  as  10  tor,  would  result  in  negligible  weight  and  pressure 
differences.  The  final  assumption  was  the  basic  premise  behind  the  RCVLTAV.  Structure 
walls  were  assumed  to  be  able  to  withstand  external  pressure  due  to  rotational  velocity 
alone. 
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™c  Pskin  '  (j-skinc  '  ®-c  ”1”  ^skins  '  &s) 


(39) 


ac  —  2  ■  n  ■  rc  ■  lc 


(40) 


as  =  4  ■  7r  ■  rc2 


(41) 


Where: 

ac  —  cylinder  area 

as  —  spherical  endcap  area 

mc  =  cylinder  mass 

The  inertial  properties  of  the  cylinder  walls  were  based  on  tskirLc  and  pskin,  which 
enabled  calculation  of  the  required  rotational  velocity  using  equation  (45).  This  equation 
was  derived  using  equations  (42)  -  (44).  The  total  force  due  to  atmospheric  pressure  (Fp) 
was  calculated  using  equation  (42)  by  multiplying  the  surface  area  of  the  cylinder  by 
Patm ■  Centripetal  force  due  to  rotational  velocity,  was  used  to  detennine  the  required 
rotational  velocity  (u>)  for  Fc  to  equal  Fp  using  equation  (43)  [43].  The  cylinder  mass  was 
calculated  using  the  surface  area  of  a  cylinder  times  the  skin  thickness  and  density  as 
shown  in  equation  (48).  The  derivation  of  equation  (45)  was  accomplished  by 
substituting  the  right  side  of  equation  (43)  for  Fp  in  equation  (42).  The  buckling  pressure 
(pcr)  of  the  stationary  cylinder  was  calculated  using  equations  (46)  -  (50),  [4,  p.  1299]. 
The  derivation  of  equations  (46)  -  (50)  are  beyond  the  scope  of  this  study,  but  can  be 
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found  in  [44].  The  amount  of  pressure  that  had  to  be  countered  using  Fc  was  reduced  by 


pcr  from  equation  (45).  It  turned  out  that  at  thickness  that  allowed  a  W/B<1,  patm  was 
about  four  orders  of  magnitude  greater  than  pcr.  Because  of  this  the  effect  of  increasing 
thickness  was  negligible  and  the  centripetal  force  Fc  was  necessary  to  counteract 
practically  the  entire  load  due  to  patm- 


Fp  2  ■  7T  ■  Vc  '  lc  Patm 


(42) 


Fc  =  mc  ■  rc  ■  a)2 


(43) 


mc  2  ■  7T  ■  rc  ‘  lc  tskinc  "  Pskir 


(44) 


u>  = 


Patm  Per 
I  Pc  '  tskinr  '  Pskir 


(45) 


Per  ®cr 


tskinc 

rc 


(46) 


KcTC  E  / tskinc\ 

~  12(1-  v2)  V  lc  ) 


(47) 
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K  =  0.6337  -  0.1455  ■  log3  z  -  0.01915  ■  log3  z 


(49) 


z  = 


l 


2 

c 


Lskinr 


(50) 


Where: 

as  —  spherical  endcap  area 

Fc  —  centripetal  force 

Fp  —  force  due  to  atmospheric  pressure 

pcr  —  critical  pressure 

or  —  radial  stress  due  to  external  pressure 

a)  —  rotational  velocity 


3.2.2  Aerodynamics 

The  cross  section  of  a  pair  of  infinitely  long  counter  rotating  cylinders  with 
streamlines  flowing  past  them  from  left  to  right  is  shown  in  Figure  22  [45,  p.  13].  The 
streamlines  on  the  left  show  cylinders  that  are  rotating  at  fl  <  flcrit,  where  f l  is  defined  in 
equation  (80)  as  the  ratio  of  cylinder  surface  speed  ur  to  freestream  velocity  um.  The 
cylinders  on  the  right  show  cylinders  rotating  at  ncrit,  which  is  the  rotational  velocity  at 
which  the  streamlines  form  a  virtual  elliptical  body  that  produces  no  vortices  and  thus 
experiences  zero  translational  drag.  Gap  size  determines  ncrit,  with  smaller  gap  sizes 
corresponding  to  smaller  flcrit.  Gap  size  is  measured  in  cylinder  diameters,  so  that  a  gap 
size  of  1  corresponds  to  two  cylinders  whose  closest  points  are  one  diameter  apart.  A  pair 


of  cylinders  with  a  gap  size  of  1  were  determined  to  both  numerically  and  experimentally 
to  have  ncrit  S3  [45,  p.  13].  The  cylinders  require  power  to  maintain  their  rotational 
velocity,  which  means  drag  is  experienced  in  the  rotational  direction,  but  not  in  the 
horizontal  direction.  In  Figure  22  the  horizontal  direction  would  be  to  the  left,  with  um 
moving  from  left  to  right.  The  streamlines  in  Figure  22  were  produced  using  a  computer 
program,  however  pictures  were  taken  during  a  physical  experiment  that  showed  the 
same  formation  of  a  virtual  elliptical  body  [45,  p.  14].  These  pictures  were  taken  at 
Re=150.  Rotational  drag  is  referred  to  as  drag  torque  (M)  because  it  is  measured  as  a 
torque  instead  of  as  a  force.  An  airship  consisting  of  a  properly  spaced  pair  of  counter 
rotating  cylinders  could  conceivably  achieve  high  forward  velocities  with  low 
translational  drag  by  forming  a  virtual  elliptical  body. 


Figure  22:  Formation  of  Virtual  Elliptical  Body  at  Critical  Velocity  |45] 


Where: 

ur  —  cylinder  tip  speed 
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itoo  =  free  stream  velocity  (translational) 


3.2.3  Power  Requirements 

The  power  required  to  rotate  a  cylinder  in  a  stationary  body  of  air  is  shown  in 
equation  (52)  and  is  the  combination  of  power  required  due  to  drag  torque  Pm  ?  spherical 
endcap  rotation  Ps,  and  bearing  friction  Pb.  Power  required  due  to  drag  torque  PM  was 
calculated  using  equation  (53)  [46],  which  is  the  rotational  velocity  times  drag  torque 
(M).  Equation  (54)  [47,  p.  370],  was  used  to  calculate  M,  where  the  rotational  drag 
coefficient  ( CDc )  was  calculated  according  to  equation  (55)  [47,  p.  370].  These  equations 
were  derived  in  [47]  based  on  experimental  analysis  of  high  speed,  rotating  cylinders  by 
Theodorson  and  aerodynamic  theory  by  Prandtl  and  von  Karman.  Interestingly, 
Theodorson  discovered  that  the  rotational  drag  coefficient,  CDc,  of  cylinders  is 
independent  of  velocity  even  though  his  experiments  ranged  well  above  cylinder  surface 
velocities  of  mach=l.The  equation  for  Reynolds  number  in  equation  (56)  [48,  p.  10]  is 
the  same  as  the  Reynolds  number  for  an  object  translating  through  air,  but  instead  of  the 
length  of  the  vehicle,  the  radius  is  used.  These  equations  were  used  to  calculate  the  power 
required  to  rotate  a  single  cylinder  in  a  stationary  body  of  air,  which  is  greater  than  the 
power  required  to  rotate  cylinders  in  a  fully  developed  virtual  elliptical  body  [45,  p.  32]. 
The  power  reduction  for  a  virtual  elliptical  body  was  considered  below,  but  it  was  first 
necessary  to  find  the  power  required  to  rotate  cylinders  without  a  fully  developed 
elliptical  body. 

Pb  was  calculated  using  equation  (63)  [49,  p.  460],  Pb  is  a  function  of  bearing 

load  ( Fb ),  bearing  radius  (rb),  <n,  and  coefficient  of  friction  (C/6).  Fb  was  assumed  to  be 

64 


the  cylinder  useful  lift,  which  is  the  buoyant  force  of  the  cylinder  minus  its  weight.  rb  and 
Cy^were  found  on  a  bearing  vendor’s  website  [50]. 

Ps  was  calculated  using  equation  (58),  [47,  p.  370]  and  was  derived  during  the 
same  study  that  produced  equations  (53)  -  (55).  The  spherical  end-caps  on  the  rotating 
cylinders  were  assumed  to  be  rotating  disks  for  drag  purposes.  The  drag  coefficient  of  a 
rotating  disk  (CmJ  was  verified  experimentally  by  [47]  and  shown  in  equation  (59). 


Pr  “  PM  +  Pb  +  ps 


(52) 


Pm  —  o  ■  M 


(53) 


M  Cqc  ■  Pair  '  Ur  ‘  TC  ■  Tc  ■  lc 


(54) 


—0.6  +  4.07  ■  logw  Re  ■  yjCDc 


(55) 


Re  Pair 


(56) 


Pb  ~  Cfb  '  rb  '  '  Fb 


(57) 
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(58) 


Ps  = 


'  Pair 


a r 


r  — 

'-'m  _ 


0.146 


Re 


0.2 


(59) 


Where: 

CDc  —  rotating  cylinder  drag  coefficient 

Cfh  —  bearing  coefficient  of  friction 

Cms  =  spherical  endcap  drag  coefficient 

Fb  —  bearing  load 

lc  —  cylinder  length 

M  —  rotating  cylinder  drag  moment 

Pb  —  power  required  due  to  bearing  friction 

PM  —  power  required  due  to  torque  drag 

Ps  —  power  required  due  to  spherical  endcap  rotational  drag 

Pr  =  total  power  required  to  rotate  cylinder 

rb  =  bearing  radius 

Re  —  Reynold's  Number 

ut  =  rotating  cylinder  tip  speed 

/r  =  air  viscosity 

a)  —  cylinder  rotational  velocity 


Figure  23,  [45]  shows  that  as  the  gap  size,  displayed  as  g  in  Figure  23,  between 
rotating  cylinders  decreases,  the  reduction  in  torque  drag  increases.  Also,  the  critical  tip 
speed  ur  required  to  achieve  the  virtual  elliptical  body  decreases  as  ncrit  decreases  as  in 
equation  (5 1).  The  implication  of  this  relationship  is  that  the  closer  the  cylinders  are 
placed  on  the  vehicle,  the  higher  translational  velocity  the  vehicle  can  achieve  while 
maintaining  a  virtual  elliptical  body  for  a  given  cylinder  rotational  velocity,  and  the  lower 
the  torque  drag  (power  required  to  rotate  cylinders)  that  will  be  encountered.  This 
research  has  not  yet  been  conducted  at  gap  sizes  of  less  than  1  cylinder  diameter.  An 


66 


interesting  area  of  future  research  would  be  to  measure  torque  and  translational  drag  on 
counter-rotating  cylinders  with  gap  sizes  between  0  and  1  diameter  to  produce  graphs 
similar  to  Figure  23. 


Re  =  200.  g*  =  1  ( Doublet  | 


— i —  Re  =  200,  g  =  3  (Doublet) 

Re  =  200,  g*  =  5  (Doublet) 

■  *  -  Re  =  200,  g'  =  1  (Reverse  doublet) 

-  +  -  Re  =  200. g’  =  3  (Reverse  doublet) 

-  ▼  -  Re  =  200,  g*  =  5  ( Reverse  doublet) 
.  Re  =  200  (Single  rotating  cylinder) 


CP  1.4  vs  1.8 
75  %  Torque  Drag 


3.2.4  Engine  and  Fuel  Weight 

After  detennining  the  translational  velocity  of  a  RCVLTAV  based  on  equation 
(5 1)  the  vehicle  payload  (W)  was  calculated  using  equation  (60)  assuming  two  counter¬ 
rotating  cylinders  with  buoyant  lift  Lcb,  a  gondola  the  size  of  a  C-27  cabin  with  mass  mg, 
4  cylinder  bearings  with  mass  mb,  a  turboshaft  engine  with  enough  power  to  rotate  the 
cylinders  and  provide  sufficient  thrust  to  overcome  the  drag  of  the  gondola  with  mass 
me,  as  well  as  enough  fuel  to  travel  2130  km  with  mass  mf.  It  was  assumed  that  no 

translational  drag  would  be  encountered  by  the  rotating  cylinders.  Total  power  required 
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(Ptot)  was  used  to  determine  fuel  and  engine  requirements.  Engine  mass  mewas 
calculated  using  equation  (61)  which  is  an  empirical  formula  for  specific  engine  power 
for  a  turbo  shaft  engine  [13,  p.  225].  Fuel  mass  mf  was  calculated  using  equation  (62) 
which  was  derived  based  on  the  work  required  to  fly  1852  km,  the  energy  density  of  fuel 
(0,  and  engine  efficiency  (kee). 

Equation  (63)  was  used  to  determine  Ptot  5  which  was  assumed  to  be  composed  of 
the  power  required  to  spin  two  cylinders  and  the  power  required  to  counter  translational 
drag  (Px).  Power  required  due  to  drag  Px  was  calculated  using  equation  (64)  which  is  the 
power  equation  for  propellers  with  thrust  equal  to  drag  [25].  This  is  a  function  of 
propeller  efficiency  (kpe),  propeller  area  (ap),  thrust  required,  which  in  this  case  is  the 
drag  due  to  the  gondola  (Dg)  assuming  steady  level  flight,  Uoo,  and  pair.  Gondola  drag  Dg 
was  calculated  based  on  gondola  cross  sectional  area  ag  (perpendicular  to  u^)  and  the 
gondola  drag  coefficient  Cdg,  [13].  The  drag  coefficient  of  the  gondola  Cdgwas  assumed 
to  be  the  same  as  the  average  subsonic  aircraft[51]. 

W  —  2  ■  Lcb  —  g  ■  (4  ■  mb  +  mg  +  me  +  mf)  (60) 

me  =  Ptot  ■  0.00201  — ^7  (61) 

watt 

PtOt 

mf  -  range  ■  - - - 

t^ee  ’  Uoo  '  S 
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(63) 


ptot  -  Px  +  2  '  Pr 


p  —  n 
rx  ^pe 


D2 

ug 


2  '  P 


air  '  '  uc 


(64) 


Un 


(65) 


Where: 

ag  —  cabin  cross  sectional  area  perpendicular  to  um 

ap  —  propeller  area 

Cd  —  gondola  drag  coefficient 

Dg  —  gondola  drag 

kee  —  engine  efficiency 

kpe  —  propeller  efficiency 

mb  =  bearing  mass 

mg  —  gondola  mass 

Ptot  =  total  power  required  for  propellers  and  rotating  cylinders 
Px  —  power  required  due  to  cabin  drag 
W  =  useful  load 
^  =  energy  density  of  fuel 


The  calculation  of  W  was  necessary  to  determine  the  feasibility  of  the  RCVLTAV 
concept.  If  W  were  negative  the  aircraft  could  not  fly.  If  W  were  positive  it  could  be  used 
along  with  um  and  Ptot  to  calculate  eVK.  Von  Karman  efficiency  could  then  be  used  to 
compare  the  performance  of  the  conceptual  vehicle  to  existing  vehicles. 
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3.3  Isogrid  Stiffened  Sphere 

An  isogrid  acts  like  an  isotropic,  thin  shelled  structure,  therefore  equation  (14)  can 
be  used  to  evaluate  an  isogrid  stiffened  sphere  using  equivalent  thickness  and  stiffness. 
The  NASA  isogrid  handbook  [36]  describes  a  simple  isogrid  composed  of  equilateral 
triangles  as  pictured  in  Figure  24,  [36,  p.  20.0.019]. 


tskins 

z? 

d 


Where: 

tskins  —  skin  thickness 
d  —  flange  height 
b  —  flange  width 
h  =  triangle  height 
a  =  triangle  leg  length 

Equation  (66)  is  the  same  as  equation  (14)  but  with  E  and  t  replaced  by  E *  and  t* 
which  are  the  equivalent  Young’s  Modulus  and  skin  thickness  respectively.  The 
equivalent  parameters  E  *  and  t*  account  for  the  stiffening  effect  of  the  grid,  but  because 
they  are  isotropic,  equation  (14)  is  still  valid.  Equations  (67)  and  (68)  give  E*  and  t* 

where  represents  stiffness  and  represents  thickness  effect  of  the  stiffeners  on 

the  skin. 


Figure  24:  Isogrid  Geometry  |36] 
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(66) 


Pcrit  0.8  ' 


E  * 

V 1  —  V2 


2 


E* 


(1  +  a)2 

E  ■- — — — 


(67) 


t  =  t 


skin? 


1  +  a 


(68) 


Where: 

E*  —  equivalent  monocoque  stiffness 
t*  =  equivalent  monocoque  thickness 

The  non-dimensional  values  B  and  a  are  given  by  equation  (69)  and  equation 
(70).  They  are  geometric  properties  of  the  isogrid  used  in  equations  (67)  and  (68).  They 
were  derived  based  on  transformed  geometry  with  the  derivation  in  [36,  p.  20.0.019]. 


P  = 


3  ■  a 


+  (1  +  a)  • 


(69) 


b  ■  d 

tskins  '  ^ 


(70) 


Equations  (66)-  (70)  can  be  used  to  determine  the  complete  sphere  buckling 
mode.  However  there  are  other  buckling  modes  that  must  be  examined.  These  are  rib  and 
skin  buckling  and  are  determined  by  equations  (71)  and  (72)  respectively  [36].  Equation 
(71)  is  the  equation  for  buckling  of  a  plate  that  is  clamped  on  three  sides  and  free  on  the 
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fourth  [36,  p.  4.1.003].  This  is  an  appropriate  model  because  the  stiffener  is  optimal  at  a/d 
<5  which  gives  it  the  dimensionality  of  a  plate.  If  a/d> 5  beam  theory  would  be  more 
appropriate.  Equation  (72)  is  the  equation  for  buckling  of  an  edge  loaded,  clamped, 
triangular  plate  [36,  p.  4.1.004]. 


n2  ■  E  (b\ 

ocrb  =  0.456  •  (1  _y2)  ■  (jJ 


5  TC2  ■  E  / tskins\ 

Gcsb~Yi\i -v*y\ir) 

Where: 

Ocrb  —  rib  buckling  critical  stress 
ocsb  —  skin  buckling  critical  stress 


(71) 


(72) 


All  values  of  t,  d,  b, and  h  that  satisfy  equations  (66),  (71)  and  (72)  are  part  of  the 
feasible  design  space  for  the  grid  stiffened  sphere.  By  simultaneously  satisfying  these 
equations,  with  the  assumption  that  the  optimal  point  would  occur  when  all  constraints  are 
active,  a  system  of  three  equations  and  four  unknowns  was  fonned.  This  allowed  a  simple 
graphical  optimization  method  in  which  the  lowest  W/B  was  chosen  by  varying  one 
variable  and  solving  for  the  remaining  three  variables.  In  this  case  tskills  was  chosen  as  the 
free  variable,  and  d,  b,  and  h  were  found  by  satisfying  equations  (66),  (71),  and  (72). 

Figure  25  is  a  plot  of  W/B  vs  tskirLs  where  W/B  is  a  function  of  smeared  skin 
thickness  (t),  rs,  and  atmospheric  properties  as  shown  in  equation  (73).  Equation  (73)  is 
the  ratio  of  sphere  weight  to  buoyant  force.  Sphere  weight  was  based  on  as,  pskin,  and  t. 
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t  was  calculated  by  smearing  the  stiffeners  and  adding  them  to  tskirLs  as  shown  in 
equation  (76).  Smearing  was  accomplished  geometrically  by  taking  stiffener  length  (a) 
times  stiffener  cross-sectional  area  ( b  ■  d),  assuming  half  of  each  stiffener  was  applied  to 
each  triangle,  assuming  three  stiffeners  per  triangle,  and  dividing  by  triangle  area  which 

was  j  ■  a  •  h.  This  gave  the  smeared  thickness  of  3  ■  b  ■  ^  for  each  stiffener. 


W 

~B 


Pskin  '  CLS  ‘  t 
Pair  '  Vs 


(73) 


t  tskins  "b  3  '  b  ■ 


(74) 


Where: 

t  —  smeared  skin  thickness 
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Figure  25:  UHM  Carbon  Epoxy  Isogrid  Optimization 


The  optimal  point  in  Figure  25  was  the  point  of  minimum  W/B.  This  W/B,  being 
a  function  of  tskiUs,  enabled  selection  of  the  optimal  variables  tski7ls,  b,  d,  and  h.  This 
method,  however,  did  not  necessarily  result  in  the  true  optimal  geometry.  It  was 
recommended  by  the  isogrid  handbook,  which  was  written  in  1973  before  more  robust 
optimization  techniques  were  made  possible  by  modern  computers.  An  alternate 
optimization  technique  was  performed  using  a  non-linear  programming  optimization 
routine  in  Matlab.  The  problem  was  fonnulated  with  equation  (73)  as  the  cost  function  to 
minimize,  equations  (66),  (71),  and  (72)  as  constraints  to  be  satisfied,  and  tskiris,  b,  d,  and 
h  as  variables.  The  optimization  program,  called  fmincon,  conducted  a  search  in  order  to 
minimize  the  cost  function  while  satisfying  the  constraints  by  systematically  changing  the 
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variables.  More  details  on  this  optimization  technique  are  available  in  3.4.2  Optimization, 
and  the  code  is  available  in  Appendix  A.  13-15. 

3.4  Geodesic  Sphere 

The  geodesic  sphere  investigated  consisted  of  a  frame  made  of  UHM  carbon  fiber 
tubes  with  a  reinforced  Mylar  membrane  stretched  over  it  similar  to  a  camping  tent. 
Figure  26  shows  an  icosahedron  where  the  black  lines  represent  the  UHM  carbon  frame, 
and  the  shaded  triangular  faces  represent  the  Vectran  reinforced  Mylar  membrane. 


Figure  26:  Geodesic  Sphere  Representation 


The  material  used  to  reinforce  the  Mylar  membrane  would  need  to  have  low 
density,  high  strength,  high  toughness  and  high  stiffness.  Two  materials  currently 
available  for  stiffening  the  skin  are  Zylon  and  Vectran  as  shown  in  Table  6.  If  graphene 
could  be  made  in  usable  sizes  and  attached  to  the  structure,  it  would  far  exceed  the 


75 


performance  of  currently  available  materials.  The  optimization  was  performed  for  the 
frame  alone.  The  skin  was  evaluated  separately  after  the  frame. 


Table  6:  Material 

Properties-Skin 

Material 

from 

Matweb 

E,  Pa 

Ps) 

kg 

m3 

Gy, 

GPa 

Geodesic 

Skin 

W/B 

Comments 

Graphene 

250 ■ 109 

100 

130 

0.02 

Only 
Available 
in  small 
sheets 

Zylon 

303  ■  109 

1560 

5.8 

2.5 

Fiber 

Vectran 

75 ■ 109 

1400 

3.2 

3.6 

Fiber 

The  W/B  of  the  UHM  carbon  tube  frame  was  computed  according  to  equation 
(75).  The  variables  tcb,  rcb,  lcb,  and  pci,  denote  the  thickness,  radius,  average  length,  and 
material  density  respectively  of  the  cylindrical  beams.  The  numerator  of  equation  (75)  is 
the  mass  of  the  frame.  It  was  derived  by  taking  the  cross  sectional  area  of  each  beam, 
multiplying  it  by  its  length  and  density  to  get  the  mass  of  each  beam,  then  multiplying  the 
mass  of  each  beam  by  the  number  of  beams  in  the  structure,  ne.  Equation  (32)  was  used 
to  calculate  ne.  The  buoyant  mass,  or  mass  of  the  air  displaced  by  the  geodesic  sphere,  is 
denoted  by  mb,  which  is  a  function  of  the  density  of  air  and  volume  of  the  geodesic 
sphere  as  shown  in  equation  (80).  An  optimization  program  called  firm  neon  was  run  in 
Matlab  to  find  the  minimum  W/B  of  the  UHM  carbon  tube  frame.  The  optimization 
problem  was  formulated  as  follows  [52]: 
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Minimize: 


F  (j-cb>  ^cb>  /) 


Pcb  '  2 ' 7T  ■  rcb  '  tcb  '  lcb  '  ?1.e 
mb 


(75) 


Subject  to: 


= 


Pe 

2mTCmY  cb'tcb 


G cbcy 


<  0 


(76) 


92 


=  Pe~ 


4  ■  n2  ■  Erh.I, 


cb'cb 


l2 

Lcb 


<  0 


(77) 


93  =  *-f  ~  rcb  <  0 


(78) 


,  ^ch 

94  =  rcb+  — 


V3 


■  n 


4  ■/ 


<  o 


(79) 


Where: 


Vgs  '  Pair 


(80) 
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And: 


g1  =  requirement  that  the  compressive  yield  strength  of  the  material  not  be 
exceeded  by  the  stress  in  any  given  beam. 
g2  =  Euler’s  buckling  equation  for  a  cylindrical  column  under  axial  load  with 
clamped  ends,  [4] 

g3  =  requirement  that  the  thickness  not  exceed  twice  the  radius,  the  thickness  of 
the  beam  cannot  extend  past  the  center  point  of  the  beam,  at  which  point  the 
cylinder  is  no  longer  a  hollow  pipe  but  is  now  a  solid  bar. 
gA  =  requirement  that  the  radius  of  the  beam  not  extend  past  the  center  of  a  given 
minor  triangle.  This  is  a  conservative  estimate  that  will  result  in  at  least  a 
small  space  between  beams  at  the  center  of  each  triangular  face. 

Icb  =  the  area  moment  of  inertia  of  the  cylindrical  beam 

Pe  =  Euler  load  on  end  of  beam 

°cbcy  —  compressive  yield  stress  of  cylindrical  beam 

Assumptions: 

Clamped  beams,  even  pressure  distribution  on  skin,  loads  on  beam  ends  only  as 
equivalent  forces,  skin  does  not  add  stiffness,  safety  factor  of  1.5,  standard  atmosphere  at 
sea  level,  failure  in  any  member  constitutes  global  failure,  no  added  mass  for  connectors 
at  ends  of  beams. 

The  constraint  gi  was  derived  based  on  the  requirement  that  the  material 
compressive  strength  of  the  cylindrical  beams,  {pcb  ),  not  be  exceeded.  This  was 
detennined  using  FEA  to  solve  for  the  stress  in  each  beam  using  analyzebeam.m  in 
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Appendix  A. 5.  After  solving  for  the  stress,  the  maximum  stress  out  of  all  the  members  in 
the  structure  was  found  and  compared  to  ocbcy.  As  long  as  <Jcbcy  was  of  greater  absolute 

value  than  the  stess  in  every  member,  gi  was  satisfied. 

Note  that  all  of  the  constraints  were  set  to  be  less  than  or  equal  to  zero.  This  was 
required  by  fmincon,  the  optimization  function  used  in  Matlab.  If  any  inequality 
constraint  was  less  than  or  equal  to  zero,  it  was  satisfied.  If  a  constraint  was  greater  than 
zero  it  was  violated  and  the  solution  was  considered  infeasible.  In  the  case  of  gi,  if 
°cbcy  was  less  than  the  stress  in  the  beam,  then  gi  would  be  positive  and  the  material 

strength  of  the  beam  would  be  exceeded. 

The  constraint  g2  was  derived  based  on  the  buckling  load  of  a  cylindrical  column 
with  clamped  ends.  Clamped  ends  were  assumed  because  all  of  the  methods  of  attaching 
the  ends  of  the  beams  to  each  other  that  were  considered  involved  a  fastening  method  that 
transmitted  bending  moments  in  addition  to  forces.  If  a  ball  and  socket  attachment  were 
made  for  the  structure,  equation  (77)  would  have  to  be  replaced  by  the  buckling  equation 
for  a  cylindrical  column  with  end  load  with  a  pinned  boundary  condition  at  each  end, 
which  is  14  the  critical  load  of  the  clamped  equation. 

The  constraint  g4  was  a  geometric  constraint  created  to  keep  the  cylindrical 
beams’  walls  from  intersecting  each  other  at  the  center  of  a  triangular  face.  It  is  based  on 
the  distance  from  the  center  of  the  triangle  to  its  edge.  Triangle  height  is  given  by 
equation  (81)  where  a  is  the  leg  length  of  a  minor  triangle.  The  distance  from  the  edge  of 
a  triangle  to  its  center  is  1/2  ■  h,  which  must  be  greater  than  or  equal  to  the  radius  plus 
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thickness  of  the  cylindrical  beams,  which  gives  g4.  Note  that  a=1.05r/f,  so  that  a  more 
conservative  solution  was  obtained  by  replacing  a  with  rs  and  dividing  by  f 


(81) 


The  area  of  the  geodesic  sphere,  ags,  was  calculated  by  summing  the  area  of  each 
triangular  face  on  the  structure,  and  although  not  used  in  the  optimization,  enabled  an 
estimate  of  skin  weight.  The  volume  of  the  geodesic  sphere,  vgs,  was  calculated  by 
summing  the  volume  of  the  triangular  pyramids  fonned  by  the  three  vertices  of  each 
triangular  face  and  the  center  point  of  the  sphere  as  illustrated  by  the  shaded  region  in 
Figure  1.  Equation  (82),  [40],  gives  the  volume  of  a  pyramid,  ( vp ),  which  is  1/3  of  base 
area,  ( ab ),  times  height  perpendicular  to  base.  In  the  case  of  a  geodesic  sphere,  the  height 
was  measured  from  the  center  of  the  sphere  to  the  centroid  of  any  triangle  on  the  surface 
and  the  base  was  the  area  of  the  triangular  face.  By  summing  the  volume  of  the  pyramids 
formed  by  every  base,  the  total  sphere  volume  was  found. 


(82) 
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Figure  27:  Geodesic  Sphere  -  Volume  Calculation 


The  Euler  Load  Pe  was  calculated  using  FEA  with  beam  elements.  Forces  on  the 
ends  of  the  beams  were  calculated  using  pair  times  the  area  of  each  triangular  face, 
distributed  over  the  three  vertices  of  the  triangle  in  the  direction  normal  to  the  face.  This 
method  resulted  in  zero  moment  on  each  beam.  This  was  because  forces  were  applied  at 
the  nodes  and  not  on  the  sides  of  the  beams.  It  was  also  because  multiple  beams  resisted 
the  equivalent  load  at  each  node,  so  that  axial  forces  in  each  beam  resisted  the  equivalent 
load  no  matter  which  way  it  was  oriented.  The  boundary  condition  consisted  of  restricting 
6  degrees  of  freedom  on  the  node  at  the  bottom  of  the  sphere.  No  mass  was  added  for 
connectors  at  the  ends  of  the  beams  based  on  the  assumption  that  the  sphere  could  be 
constructed  of  continuously  wound  fdaments  about  a  mandrill  where  vertices  would  be  a 
flat,  continuous  transition  between  pultruded  rods.  If  connectors  are  required  in  the  future 
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their  mass  will  have  to  be  added  to  the  cost  function  and  would  be  the  weight  of  a 
connector  multiplied  by  the  number  of  vertices,  nv. 

Table  7  displays  the  material  properties  used  to  calculate  the  optimal  sphere 
geometry.  Additional  constants  used  include  rs  =  0.33  m,  pa=  101325  Pa,  and  pa  — 

k  g 

1.225  — .  The  assumption  was  made  that  UHM  carbon  had  a  1/^=0. 6  of  the  composite 

tubes  and  the  composite  stiffness  in  the  longitudinal  direction  was  calculated  using 
equation  (24).  Poisson’s  ratio  was  assumed  to  be  0.3. 


Table  7:  Material  Properties-Frame 


Material 

E,  Pa 

„  kg 

P  ’  m3 

V 

®CbCy>  PCL 

Dupont  E-130-X  Pitch 
Based  UHM  Carbon 
Fiber 

895 ■ 109 

2130 

Unpublished 

Unpublished 

UHM  Unidirectional 

Carbon  Epoxy  tubes 

538-  109 

1522 

0.30 

1  ■  109 

Figure  28  shows  a  workflow  diagram  of  the  optimization  process.  All  programs 
were  run  within  Matlab  to  include  the  geometry,  force  calculations,  FEA,  and 
optimization  routine.  Optimization  was  conducted  graphically  as  well  as  numerically. 
The  graphical  method  was  primarily  used  to  validate  the  numerical  results,  as  well  as 
understand  the  design  space.  The  numerical  method  was  conducted  in  LTAopt.m, 
Appendix  A.  1 ,  which  used  the  Matlab  optimization  routine  fmincon.  fmincon  is  a 
constrained,  nonlinear  program  that  worked  by  systematically  changing  the  variables  f 
rcb,  and  tch  in  order  to  minimize  the  cost,  which  was  evaluated  using  Cost.m,  Appendix 
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A.2,  where  cost  was  defined  by  equation  (75).  Cost. m  called  the  Matlab  program 

Sphere_param.m,  Appendix  A. 3,  which  used  the  variables  from  fmincon  to  calculate 

sphere  parameters  that  were  then  used  by  Cost.m  to  calculate  the  W/B  of  the  frame.  After 

detennining  that  it  had  moved  in  a  direction  to  minimize  the  cost,  fmincon  checked  the 

feasibility  of  the  new  variables,  fmincon  can  automatically  check  the  feasibility  of  linear 

constraints,  which  it  did  for  g3,  but  it  must  call  a  function  to  evaluate  nonlinear 

constraints.  It  did  this  by  calling  Constraints.m,  Appendix  A. 4.  Constraints. m  was  used  to 

calculate  gi,  g2,  and  g4.  g4  was  simply  a  nonlinear  function  of  the  input  variables,  and  so 

was  quickly  evaluated  by  Constraints.m,  but  gl  and  g2  required  the  stress  in  each 

member  to  be  evaluated,  which  required  geometry  and  FEA  based  on  the  input  variables 

from  fmincon.  Constraints.m  called  Analyzebeam.m  [53],  Appendix  A. 5,  in  order  to 

retrieve  the  stress  in  the  members.  Analyze_beam.m  detennined  the  stresses  in  the 

members  by  retrieving  a  finite  element  input  card,  sending  it  to  a  finite  element  reader 

and  evaluating  the  resulting  finite  element  model  using  FEA.  Analyze  beam.m  retrieved 

the  input  card  by  calling  the  Matlab  program  Datcreate.m  [54],  Appendix  A. 6,  which 

created  a  .dat  input  file.  Datcreate.m  created  the  input  file  using  a  node  matrix,  element 

matrix,  force  matrix,  and  boundary  conditions  provided  by  the  Matlab  program 

LTAgeod2.m,  Appendix  A. 7.  LTAgeod2.m  calculated  the  node  matrix,  element  matrix, 

and  force  matrix  based  on / using  geodesic  geometry.  This  program  was  also  used  to 

generate  the  geodesic  figures  shown  in  this  document.  It  relied  on  the  geometry  of  an 

icosahedron  which  was  discussed  in  detail  in  chapter  2.5  Geodesic  Sphere,  as  well  as 

geodesic  geometry  discussed  later  in  this  chapter.  After  Datcreate.m  retrieved  the 

geometric  data  from  LTAgeod2.m,  it  compiled  the  data  into  a  .dat  input  file  for  use  by 
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Analyzebeam.m.  After  receiving  the  .dat  input  file,  Analyzebeam.m  called  the  .dat  file 
reader,  Datread_beam.m  [55],  Appendix  A. 8,  which  processed  the  .dat  file. 

Analyze  beam.m  then  created  a  global  stiffness  matrix,  force  vector,  displacement 
vector,  and  solved  for  displacements.  This  data  was  then  used  to  determine  stress  in  each 
member  of  the  frame,  which  was  used  by  Constraints. m  to  calculate  gi  and  go .  Finally, 
fmincon  evaluated  the  constraints  to  ensure  all  were  satisfied.  It  also  checked  to  see  if  the 
change  in  cost  was  low  enough  to  indicate  an  optimal  point  had  been  found.  If 
optimization  was  not  achieved  fmincon  would  change  the  input  variables  using  the 
interior  -point  search  method  and  the  whole  process  was  repeated. 


Figure  28:  Workflow  Diagram 
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3.4.1  Geometry 

The  geodesic  sphere  geometry  was  calculated  using  the  LTAgeod2.m  program, 

Appendix  A. 7.  This  program  generated  a  geodesic  sphere  by  subdividing  an  icosahedron 

along  its  major  triangles  by  the  integer/ according  to  the  geometry  described  in  [37]. 
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Geodesix  sphere  geometry  was  calculated  by  systematically  working  from  top  down  on 
the  icosahedron.  Each  of  the  major  triangles  were  divided  into  three  groups;  the  top 
pentagon  cap  of  five  triangles  that  touch  vertex  1  in  Table  5,  the  center  band  of  10 
triangles,  and  the  bottom  pentagon  cap  of  five  triangles  touching  vertex  12  in  Table  5. 
The  major  triangles  within  these  groups  were  subdivided  by  dividing  each  major  triangle 
vertically/ times  per  major  triangle.  Planer  rings  were  then  considered  working  from  the 
top  vertex  down,  where  each  ring  consisted  of  every  vertex  on  a  plane  of  constant  <p 
(parallel  to  the  xy  cartesian  plane)  as  shown  in  Figure  29.  For  the  top  pentagon  cap  each 
subsequent  ring  had  one  more  vertex  than  the  one  before  it  for  each  major  triangle.  This 
pattern  enabled  automatic  generation  of  vertex  points  for  the  top  pentagon  cap  for  any  / 
This  pattern  is  shown  in  Figure  29,  where  each  number  represents  the  order  of  vertex 
generation  for  a  given  planer  ring  for  a  geodesic  sphere  of f=5.  Each  vertex  on  a  plane 
had  a  constant  <p.  9  was  detennined  by  dividing  2  ■  7r/5  by  the  number  of  divisions  per 
major  triangle  for  a  given  planer  ring.  As  shown  in  Figure  29,  the  number  of  divisions 
increased  by  one  for  each  planer  ring,  with  the  bottom  ring  having/ divisions.  Note  that 
the  major  triangle  shown  in  Figure  29  is  divided  into  5  rows  of  triangles  corresponding  to 
/=  5.  Also  note  that  there  are  six  triangle  edges  emenating  from  each  vertex,  except  for 
those  vertices  at  the  comers  of  major  triangles,  these  have  five  triangle  edges  and  thuse 
make  pentagon  caps. 
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Figure  29:  Geodesic  Sphere  Generation-Top  Pentagon  Cap 


The  center  band  of  10  triangles,  as  shown  in  Figure  30,  continued  the  pattern  of  <fi 
generation  for  each  planer  ring  by  dividing  the  triangles  vertically / times,  just  as  was 
done  for  the  top  pentagon  cap.  9  was  different,  however,  in  that  every  planer  ring  had 
exactly  the  same  number  of  divisions  for  the  center  band  of  10  triangles.  This  occurred 
because  the  10  triangles  along  the  center  band  alternate  between  point  up  and  point 
down.  As  one  triangle  gets  wider,  the  adjacent  triangles  get  narrower.  If  layed  flat,  the 
center  band  of  triangles  would  make  a  parallelagram  with  sides  of  length  —  5  ■  a  on  top 
and  bottom  and  length=a  on  the  sides,  where  a  is  the  length  of  a  major  triangle  leg. 
Equally  dividing  this  parallelagram  for  each  planer  ring  resulted  in  increments  of 
8  —  (2  ■  7t)/(5  ■  /).  This  pattern  again  allowed  for  automatic  generation  of  vertices  based 
on  f. 
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Figure  30:  Geodesic  Sphere  Generation-Center  Band  of  10  Major  Triangles 


The  bottom  pentagon  cap  was  simply  the  reverse  of  the  top  pentagon  cap,  which 
again  allowed  for  automatic  vertex  generation.  Similar,  although  somewhat  more 
complicated  approaches,  were  used  to  determine  beam  connections  between  vertices. 

Forces  were  calculated  by  multiplying  the  area  of  each  triangle  by  the  pressure  of 
atmosphere  and  a  1.5  safety  factor.  The  resulting  force  was  divided  by  3  and  applied  at 
the  vertices  of  each  triangle.  For  the  icosahedron,  which  consists  entirely  of  equilateral 
triangles,  this  resulted  in  completely  symmetric  loading.  For  geodesic  spheres  of />  1  this 
resulted  in  uneven  loading  due  to  slightly  different  triangle  areas.  This  caused  some 
beams  to  have  higher  stress  than  others. 

3.4.2  Optimization 

The  fmincon  optimization  routine  in  Matlab  is  a  robust  solver  that  finds  the 
minimum  value  of  a  constrained,  nonlinear,  multivariable  function.  For  an  in  depth 
explanation  of  the  theory  behind  fmincon  see  the  Matlab  help  file,  which  includes 
references  that  explain  the  theory  [56].  In  this  case,  equation  (75)  was  minimized  subject 
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to  the  constraints  gl,  g2,  g3,  and  g4.  The  constraints  gl  and  g2  required  outputs  from  the 
FEA  program.  The  FEA  program  was  limited  in  function  to  6  degree  of  freedom  beam 
elements,  but  it  was  computationally  efficient  and  operated  within  Matlab,  which 
significantly  increased  the  speed  of  each  iteration  above  what  would  have  been  possible 
using  commercial  FEA  codes.  The  results  of  the  FEA  program  were  checked  against 
NASTRAN,  which  by  design  uses  the  same  input  files  as  the  FEA  program  used  for  this 
study.  Stress  results  between  NASTRAN  and  Analyze_beam.m  were  the  same  when 
given  the  same  .dat  file  from  Datcreate.m. 

A  genetic  algorithm  minimization  routine  was  also  used  in  Matlab,  called  ga.m,  to 

make  sure  fmincon  was  not  being  caught  in  local  minimums.  The  genetic  algorithm 

operates  by  testing  cost  and  feasibility  at  multiple  points  in  the  design  space  over 

successive  generations  of  design  point  populations.  The  first  generation  is  often  random, 

but  subsequent  generations  are  based  on  the  best  design  points,  tenned  individuals,  from 

the  previous  generation.  After  testing  these  points,  and  assigning  a  fitness  value  to  them, 

the  genetic  algorithm  selects  the  most  fit  individuals,  and  reproduces  them.  The 

subsequent  generation  is  formed  by  a  random  group  of  individuals  that  are  based  on  the 

most  fit  individuals  from  the  previous  generation,  similar  to  the  theory  of  natural 

selection.This  continues  until  the  genetic  algorithm  converges  sufficiently  or  reaches  the 

specified  number  of  generations.  This  method  was  highly  inefficient  because  it  required 

multiple  iterations  of  the  optimization  process  shown  in  Figure  28  for  each  generation, 

and  it  required  multiple  generations.  The  purpose  of  using  the  genetic  algorithm  was  that 

it  was  less  susceptible  to  getting  stuck  in  a  local  minimum  due  to  use  of  the  wrong  start 

point.  Unlike  fmincon,  which  requires  an  initial  value  for  each  variable,  the  genetic 
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algorithm  only  required  upper  and  lower  bounds  on  variables.  The  optimization  results  of 
ga.m  were  nearly  identical  to  fmincon.  fmincon  was  much  faster,  and  more  consistent 
however.  The  consistency  of  fmincon  can  be  explained  by  its  search  method,  which  uses 
first  and  second  derivative  infonnation  .  The  genetic  algorithm  used  a  more  random 
approach,  in  which  no  two  searches  tested  the  exact  same  populations  of  points,  so  that 
the  optimal  individual  was  always  a  little  different.  This  imprecise  nature  of  genetic 
algorithms  was  acceptable  as  long  as  the  number  of  generations,  and  number  of 
individuals  per  generation  were  sufficient.  The  tolerances  on  fitness  value  change  over 
subsequent  generations  could  be  adjusted  so  that  the  genetic  algorithm  would  keep 
searching  until  it  began  finding  solutions  within  the  specified  error.  This  resulted  in 
satisfactory  results  within  tolerance,  but  the  solution  was  still  slightly  different  every  time 
the  program  was  run. 

Finally,  a  graphical  method  of  optimization  was  perfonned  where  the  cost 
function  and  constraint  functions  were  plotted  for  various  cylinder  radii  and  thicknesses 
while  holding/ constant.  This  provided  a  good  sanity  check  for  the  answers  fmincon 
produced  as  well  as  a  better  understanding  of  the  design  space.  An  example  of  the  design 
space  fonned  using  the  graphical  method  is  shown  in  Figure  43.  The  cost  was  a  function 
of  three  variables;  f,  rcb,  and  tcb  which  would  require  a  four  dimensional  plot  to  optimize 
graphically.  Four  dimensions  is  difficult  if  not  impossible  to  plot  or  even  comprehend.  A 
plot  of  cost  as  a  function  of  two  variables,  however,  would  only  require  three  dimensions. 
Due  to  the  long  times  required  to  generate  geometry  and  evaluate  stresses  using  FEA  for 
high //of  only  1  through  10  were  evaluated.  It  was  hoped  that  by  evaluating  ten 


geometric  frequencies  a  pattern  would  emerge,  either  high/ would  produce  a  lower  W/B 
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or  a  higher  W/B.  Because  the  variable / consisted  of  integers,  it  was  possible  to  produce 
ten  three  dimensional  plots  of  the  design  space,  each  plot  having  a  constant  f  with  cost  as 
a  function  of  rcb  and  tcb  only.  These  plots  allowed  visual  understanding  of  the  design 
space  and  provided  validation  of  fmincon  and  genetic  algorithm  results.  These  plots  can 
be  seen  in  Figure  43  and  Figure  41  and  are  explained  in  chapter  4. 

3.5  Summary 

The  methods  used  in  this  investigation  included  buoyancy,  weight,  aerodynamic, 
and  power  requirements  for  the  rotating  cylinder  vehicle,  isogrid  calculations  for  the 
isogrid  stiffened  structure,  and  geometry,  FEA  and  optimization  for  the  geodesic 
structure. 


IV.  Analysis  and  Results 

4.1  Chapter  Overview 

This  chapter  gives  the  results  of  each  study  described  in  the  methodology  section. 
Also  introduced  are  results  to  minor  investigations  which  were  necessary  throughout  the 
course  of  this  investigation  to  focus  the  research.  These  include  the  identification  of 
important  material  properties  and  choosing  an  altitude  for  evaluating  vacuum  LTA 
vehicles.  The  W/Bs  obtained  for  the  grid  stiffened  and  geodesic  spheres  are  listed  and 
were  evaluated  for  their  impact  on  potential  missions.  Plots  of  von  Karman  efficiency 
were  used  to  give  a  visual  representation  of  these  results. 
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4.2  Minor  Investigations 

The  main  focus  of  this  investigation  was  development  of  a  structure  that  could 
achieve  positive  buoyancy  using  a  vacuum.  Whenever  an  idea  for  a  structure  was 
evaluated  it  became  necessary  to  determine  the  ideal  material  for  that  structure  in  order  to 
find  the  lowest  W/B.  Reduction  of  W/B  required  improved  geometric  stiffening  or 
improved  material  properties.  The  effect  of  increasing  isotropic  stiffness  or  decreasing 
density  while  keeping  all  other  material  properties  constant  is  displayed  in  Figure  3 1 .  The 
dominant  material  property  is  E /p2,  therefore  reducing  density  has  a  greater  effect  than 
increasing  stiffness.  After  an  extensive  material  search  it  became  apparent  that  most 
materials  with  low  densities  and  high  stiffness  also  have  low  compressive  strength,  so 
that  compressive  strength  becomes  the  dominant  material  property.  These  materials  are 
generally  foams  or  micro-lattice  structures.  Carbon  epoxy  composites  appear  to  be  the 
most  promising  materials  available  for  this  problem  due  to  their  combination  of  high 
specific  stiffness  and  high  compressive  strength.  There  may  be  unpublished  materials 
with  higher  specific  stiffness  than  evaluated  in  this  study.  The  question  was  asked,  what 
happens  if  there  are  better  materials,  and  what  would  the  effect  of  improved  material 
properties  be  on  W/B?  Figure  3 1  shows  the  linear  relationship  between  density  and  W/B. 
The  assumptions  made  for  creating  this  graph  were  an  isotropic  stiffness  of  E=895  GPa, 
v  =  0.3,  oCy  =  1.7  GPa  and  using  an  isogrid  stiffened  sphere.  No  material  on  the  market 

has  these  characteristics,  but  they  serve  to  show  the  effect  of  changing  material  properties 
on  achievable  W/B  in  an  isogrid  stiffened  sphere.  These  effects  are  specific  to  the  isogrid 
geometry,  however  similar  affects  will  be  seen  on  the  frame  of  the  geodesic  sphere. 
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Figure  31:  Effect  of  Density  on  W/B 


Figure  32  shows  the  effect  of  increasing  stiffness.  Although  increasing  stiffness 
decreases  W/B,  the  W/B  is  less  sensitive  to  stiffness  than  it  is  to  density.  This  makes 
sense  when  looking  at  equation  (23),  which  shows  that  W/B  is  a  function  of  E / p2 .  The 
assumptions  made  for  creating  this  graph  were  p  —  1522  kg /m3,v  —  0.3,  and  using  an 
isogrid  stiffened  sphere.  Figure  3 1  and  Figure  32  were  produced  using  LTAsphere.m, 
Appendix  A.  10 
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Figure  32:  Effect  of  Young’s  Modulus  of  Elasticity  on  W/B 
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A  study  was  conducted  to  determine  the  optimal  internal  pressure  for  air  in  a 
partial  vacuum  LTA  structure.  The  question  that  drove  this  study  was  whether  there  was 
an  optimal  internal  air  pressure  at  which  the  required  structure  weight  due  to  external 
pressure  compared  most  favorably  with  a  decrease  of  buoyant  force  due  to  internal  air 
pressure.  Figure  33  shows  a  plot  of  W/B  at  sea  level  for  a  blade  stiffened  sphere 
assuming  it  is  designed  for  a  specific  internal  pressure  corresponding  to  the  value  on  the 
abscissa.  This  means  that  as  internal  pressure  increases  on  the  x  axis,  the  structure  is 
designed  to  withstand  less  external  pressure  loading.  As  design  internal  pressure 
increases,  even  with  the  reduced  stiffness  requirement,  the  increase  in  weight  due  to  the 
air  molecules  inside  the  structure,  as  well  as  the  decrease  in  buoyant  force  results  in  a 
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decrease  in  achievable  W/B.  These  calculations  were  perfonned  in  the  Matlab  program 
titled  LTApvalt.m  ,  Appendix  A.l  1.  The  governing  equations  in  LTApvalt  are  equations 
(66)  through  (70)  where  t,  t*,  and  E*  were  calculated  in  the  Matlab  program  titled 
LTAsphere.m  for  a  blade  stiffened  sphere. 

Figure  33  was  generated  assuming  an  isogrid  structure  with  an  isotropic  material 
with  E=450  GPa,  p  =  1522  kg/m3,  v  =  0.3,  sea  level  standard  day  conditions. 

Although  a  material  with  these  properties  does  not  necessarily  exist,  Figure  33  does  show 
the  appropriate  relationship  between  internal  pressure  and  W/B.  Note  that  the  curve  in 
Figure  33  appears  to  show  a  decrease  in  W/B  at  an  internal  pressure  of  about  80%  of  sea 
level  pressure.  There  is  indeed  a  decrease  in  W/B,  but  this  portion  of  the  curve  will 
always  decrease  to  a  minimum  W/B=l  at  an  internal  pressure  of  sea  level.  This  is  because 
the  required  structural  stiffness  approaches  zero  as  internal  pressure  approaches  sea  level 
pressure,  but  the  weight  of  air  inside  the  structure  approaches  the  weight  of  air  displaced, 
so  that  as  the  structure  weight  decreases,  the  internal  weight  of  air  increases.  The  farthest 
right  point  on  the  graph  is  where  W/B=l  because  the  structure  no  longer  exists,  and  the 
air  is  neutrally  buoyant  within  itself.  Therefore  the  only  truly  informative  part  of  the 
curve  in  Figure  33  is  the  portion  of  the  graph  to  the  left  of  the  W/B  maximum  point.  This 
portion  of  the  curve  shows  that  as  internal  pressure  decreases,  W/B  also  decreases,  so  that 
the  point  of  minimum  W/B  will  be  achieved  when  there  is  a  perfect  vacuum  inside  the 
structure. 
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Figure  33:  Effect  of  Increasing  Internal  Pressure 


Does  an  increase  in  volume  result  in  a  decrease  in  W/B?  This  question  is  often 
accompanied  with  a  reference  to  the  squared  cube  law.  The  squared  cube  law  states  that 
as  the  radius  of  a  sphere  increases,  its  surface  area  grows  at  a  rate  of  radius  squared, 
while  its  volume  grows  at  a  rate  of  radius  cubed.  These  relationships  can  be  seen  by 
comparing  equations  (21)  and  (41),  which  are  the  volume  of  a  sphere  and  surface  area  of 
a  sphere  respectively.  This  is  also  reflected  by  von  Karman  efficiency,  defined  by 
equation  (2),  which  favors  an  increase  in  volume  for  an  airship.  This  is  shown  graphically 
in  Figure  9  where  the  airship  with  the  larger  radius  (and  thus  volume)  has  a  better  von 
Karman  efficiency  (its  curve  is  down  and  to  the  right).  This  is  because  the  lifting  capacity 
of  the  airship  grows  faster  with  an  increase  in  volume  than  the  drag  of  the  airship. 
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If  the  square  cubed  law  makes  the  conventional  airship  favor  larger  sizes,  does  it 
do  the  same  to  the  vacuum  structure?  As  far  as  von  Kannan  efficiency  is  concerned,  the 
vacuum  airship  would  perfonn  better  at  higher  volumes.  But  as  far  as  W/B  is  concerned, 
the  squared  cube  law  does  not  apply.  This  is  because  the  vacuum  LTA  vehicle  must  be 
able  to  resist  compressive  loads  without  buckling  or  exceeding  material  strength.  This 
can  only  be  accomplished  by  an  increase  in  thickness  whenever  there  is  an  increase  in 
volume  for  a  given  design.  The  increase  in  thickness,  accompanied  by  the  increased 
surface  area,  increases  the  volume  of  the  skin  material  at  the  same  rate  as  the  volume  of 
the  overall  structure  is  increased,  so  that  the  mass  of  the  structure  maintains  the  same 
ratio  with  the  mass  of  the  air  displaced.  In  other  words,  the  W/B  is  constant  for  a  given 
design  regardless  of  volume.  There  are  several  advantages  to  making  a  vacuum  LTA 
vehicle  larger,  however.  One  is  that  a  larger  structure  has  thicker  walls.  The  designs  in 
this  study  require  very  thin  walls,  some  that  might  be  thinner  than  current  manufacturing 
capabilities  allow.  By  making  a  larger  structure  it  would  be  possible  to  use  larger 
thicknesses  that  can  be  accommodated  by  current  technology.  The  second  advantage  is 
increased  von  Kannan  efficiency,  as  previously  discussed.  The  property  of  constant  W/B 
regardless  of  sphere  size  means  that  the  dimensions  of  the  results  of  this  study  can  be 
scaled,  but  the  optimal  W/B  will  be  the  same  for  any  sized  sphere. 

4.3  Rotating  Cylinder  Vacuum  LTA  Vehicle 

The  RCVLTAV  was  evaluated  with  assumptions  listed  in  Table  8.  The  results 
gave  energy  requirements  so  high  for  each  radius  and  length  calculation  that  the  weight 
of  the  engine  and  fuel  was  prohibitive  for  flight.  This  was  because  the  entire  lift  capacity 
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of  the  vehicle  was  less  than  the  weight  of  the  engine  and  fuel  required  to  rotate  the 
cylinders  according  to  equations  (61)  and  (62).  This  held  true  for  rotating  cylinder 
vehicles  at  all  altitudes  with  cylinder  radii  of  10  meters  or  less  and  a  length/diameter 
<500  even  if  the  weight  of  fuel  was  neglected.  The  power  requirements  of  the 
RCVLTAV  are  plotted  with  the  power  requirements  of  the  C-27J  in  Figure  34  to  give  an 
idea  of  how  much  power  the  RCVLTAV  would  require.  An  altitude  of  19.8km  (65k  ft) 
MSL  was  assumed,  which  was  the  most  efficient  altitude  for  the  RCVLTAV.  These  high 
power  requirements  were  due  to  high  rotational  velocities,  which  were  required  by  the 
assumption  that  the  external  pressure  due  to  atmosphere  was  fully  resisted  by  the  inertia 
of  rotating  cylinder  walls.  At  sea  level  the  cylinder  surface  velocity,  ur  =  730m/s, 
which  is  greater  than  the  speed  of  sound,  while  at  19.8km  MSL  it  was  about  195  m/s. 
Due  to  the  fact  that  CDc  does  not  change  with  Mach  number  [47],  equation  (55)  remains 
valid  even  for  high  velocities.  Whether  or  not  counter  rotating  cylinders  would  form 
virtual  elliptical  bodies  at  high  Reynolds  numbers  is  uncertain.  The  RCVLTAV  at 
19.8km  would  have  Reynolds  numbers  ranging  from  1  ■  107to  1.2  ■  108  for  radii  of  1  to 
10  meters.  The  assumption  was  made  that  the  rotating  cylinders  would  act  in  the  same 
manner  as  rotating  cylinders  with  low  Reynolds  numbers.  If  the  results  of  this  study  was 
favorable  to  the  performance  of  the  RCVLTAV  this  assumption  would  need  to  be 
verified. 

Table  8:  Rotating  Cylinder  Vehicle  Assumptions 

a,g  —  10.5  m2  based  on  cylindrical  cross  section  of  C-27J  cabin[6] 
ap  —  12.6  m2  based  on  historical  airship  rotor  sizes  [9] 

Cdg  —  0.012  based  on  wind  tunnel  testing[13] 
kee  =  0.33  based  on  turbine  engine  efficiency  [13] 
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kpe  —  0.55  assuming  stern  mounted  propellers  [24] 
mb  —  26  kg  from  vendor  website  [50] 

mg  —  2114  kg  based  on  1 1  ^  and  volume  of  C-27J  cabin  [13], [6] 
range=  1852  km  based  on  C-27J  range  [6] 
l  =  43.02-  106- BP  Jet  A  [57] 

kg 

n  =  1.5  assuming  less  than  1  cylinder  diameter  spacing  [45] 
sea  level  to  19.8km  (65k  ft)  MSL  standard  day  atmospheric  conditions 


Figure  34:  Power  Required  by  RCVLTAV  at  75%  Drag  Torque  at  19.8km  MSL 


A  reduction  in  drag  torque  of  about  25%  on  counter  rotating  cylinders  due  to  a 
virtual  elliptical  body  has  been  shown  [45].  There  is  an  apparent  trend  that  as  the  gap 
between  cylinders  is  decreased,  the  drag  torque  also  decreases.  The  question  was  raised, 
how  much  would  the  drag  torque  have  to  decrease  in  order  for  the  RCVLTAV  to  become 
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feasible?  The  energy  requirements  of  a  RCVLTAV  with  a  99%  decrease  in  drag  torque, 
at  the  optimal  operating  altitude  of  19.8km  MSL,  resulted  in  a  vehicle  that  could  fly.  The 
energy  comparison  between  this  vehicle  at  various  cylinder  radii  and  the  C-27J  is  plotted 
in  Figure  35.  This  shows  that  with  a  radius  of  less  than  nine  meters,  the  RCVLTAV 
would  use  less  energy  than  the  C-27J.  In  order  to  better  understand  the  transportation 
efficiency  of  the  RCVLTAV  a  plot  of  von  Kannan  Efficiency  is  shown  for  radii  of  one 
through  ten  meters  in  Figure  36. 


Figure  35:  Power  Required  by  RCVLTAV  at  1%  Drag  Torque  at  19.8km  MSL 
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Figure  36:  RCVLTAV  at  1%  Drag  Torque  at  19.8km  MSL 


The  von  Kannan  efficiency  of  the  RCVLTAVs  in  Figure  36  are  excellent,  but  it  is 
unlikely  that  a  drag  torque  reduction  of  99%  will  occur  for  closely  spaced,  or  even 
touching  cylinders.  Even  if  a  drag  torque  reduction  of  99%  were  possible,  the  vehicle  still 
could  not  operate  at  sea  level.  At  sea  level,  the  energy  requirements  were  higher  due  to 
increased  drag,  so  that  an  engine  that  could  provide  the  required  power  would  be  too 
heavy.  Therefore  the  vehicle  would  have  to  be  assisted  in  its  climb  to  altitude  by  another 
vehicle  capable  of  carrying  itself  and  the  vehicle  in  question.  The  geometry  of  this 
structure  would  be  1  to  10  meters  in  radius  (or  higher  if  desired),  W/B=0.5 1,  305  to  3100 
meters  long,  and  have  a  mass  of  5,865  to  990,000  kg  (depending  on  radius). The 
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dimensions  of  this  vehicle,  the  inability  to  fly  at  sea  level  atmospheric  conditions,  and  the 
requirement  for  a  99%  reduction  in  drag  torque,  all  combine  to  make  the  RCVLTAV 
appear  infeasible. 

4.4  Isogrid  Stiffened  Sphere 

Using  equations  (66)  -  (72)  the  optimal  isogrid  dimensions  were  calculated  for  a 
1/3  meter  radius  sphere  made  of  the  UHM  carbon  epoxy  with  properties  shown  in  Table 
4.  The  W/B  of  various  structures  were  found  to  be  independent  of  size  so  that  any  radius 
could  have  been  chosen.  A  1/3  meter  radius  was  chosen  to  create  a  structure  capable  of 
fitting  through  a  door.  This  had  both  research  and  mission  purposes.  For  the  purposes  of 
research  it  would  be  more  economical  to  construct  a  structure  that  uses  less  material.  If  a 
vacuum  LTA  structure  could  be  designed,  it  would  be  advantageous  to  have  it  small 
enough  to  move  in  and  out  of  labs,  class  rooms,  vehicles  etc.  Also,  one  mission  type 
identified  as  suitable  for  the  vacuum  LTA  vehicle  was  as  a  small  sensor  platform  for  use 
in  the  urban  environment. 

An  optimization  was  performed  using  equations  (66),  (71),  and  (72)  and  four 
unknowns  ( b ,  d,  h,  and  tskinJ.  The  optimization  was  performed  by  assuming  the 
optimum  design  occurred  when  each  of  the  constraints  were  simultaneously  satisfied. 
Because  there  was  one  more  variable  than  equation,  it  was  possible  to  solve  the  equations 
for  multiple  values  of  tskiTls .  The  equivalent  mass  thickness,  t,  was  calculated  using 
equation  (74),  which  gave  the  thickness  that  would  be  obtained  if  the  isogrid  was 
smeared  evenly  over  the  skin. 
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Equations  (20),  (21),  (73),  and  (74)  were  used  to  determine  the  W/B  of  the  grid 
stiffened  sphere  at  standard  day  sea  level  pressure  with  a  safety  factor  of  1.5.  By  varying 
thickness  and  solving  for  the  remaining  three  variables  a  different  W/B  was  calculated  for 
each  thickness  and  plotted  in  Figure  37  for  a  sphere  made  of  Beryllium.  The  thickness 
corresponding  to  the  lowest  W/B  was  chosen  as  the  optimal  structure.  The  W/B  for  the 
optimal  structure  came  out  to  be  1.7.  Unfortunately  this  isogrid  structure  would  not 
achieve  positive  buoyancy. 


Figure  37:  Beryllium  Isogrid  Optimization 
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Equations  (5)  and  (21)  were  then  used  to  determine  the  W/B  ratio  of  an 
unstiffened  sphere  of  the  same  radius,  material,  and  atmospheric  conditions.  This  ratio 
came  out  to  be  3.6.  Therefore  the  isogrid  structure  reduced  the  W/B  ratio  of  the 
unstiffened  sphere  by  about  52%.  Both  values  can  be  seen  on  Figure  37.  This  exercise 
demonstrates  the  difficulty  of  achieving  positive  buoyancy  with  a  vacuum  structure.  It 
also  demonstrates  the  importance  of  both  material  selection  and  geometry.  Beryllium  was 
selected  because  it  had  the  highest  E/p  2  value  of  any  isotropic  material  investigated.  Note 
the  flat  portion  of  the  curve  at  the  bottom  of  Figure  37  where  W/B  is  1.7  for  a  range  of 
thicknesses.  The  flat  portion  is  caused  by  the  material  compressive  strength  constraint 
becoming  active  as  opposed  to  the  buckling  constraint.  When  this  occurs  there  is  a 
minimum  cross  sectional  area  required  in  order  to  prevent  material  failure,  so  no  matter 
what  combination  of  variables  are  chosen,  they  must  add  up  to  a  minimum  t.  If  the 
compressive  yield  strength  of  the  beryllium  could  be  augmented,  a  lower  value  for  W/B 
could  be  achieved.  Augmentation  in  this  case  could  take  the  fonn  of  adding  another 
material  to  the  structure.  For  a  composite  material  it  could  mean  adding  a  different  fiber 
or  additive  to  the  matrix.  This  demonstrates  that  although  E/p2  is  important,  in  this  case 
the  material  compressive  yield  strength  is  the  active  constraint  at  the  optimal  point,  so 
that  an  increase  in  material  stiffness  would  have  no  effect  on  W/B.  There  are  lightweight 
materials  such  as  plastic  foams  that  have  better  E/p2  values  than  UHM  carbon  fiber,  but 
due  to  their  low  compressive  yield  strengths  they  are  unsuitable  for  this  problem.  These 
materials  fail  when  their  compressive  yield  strength  is  exceeded. 

An  evaluation  of  the  grid  stiffened  sphere  was  conducted  using  Dupont  E-130-X 


carbon  fiber  in  an  Epoxy  matrix.  A  unidirectional  arrangement  of  fibers  was  assumed  for 
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the  stiffeners,  and  a  quasi-isotropic  laminate  for  the  skin.  If  a  [0/+45/90]s  ply 
arrangement  were  used  for  the  skin,  the  estimated  isotropic  stiffness  would  be 
approximately  202  GPa.  The  unidirectional  fiber  Modulus  in  the  stiffeners  was 
approximately  538  GPa.  In  order  to  use  the  isotropic  assumptions  of  the  isogrid 
handbook,  a  transformed  thickness  had  to  be  computed  for  the  stiffener.  This  was  done 
using  equation  (83),  which  gives  equivalent  stiffener  thickness  by  multiplying  the  actual 
thickness  by  the  ratio  of  stiffener  stiffness  to  skin  stiffness.  The  transformed  thickness  b* 
is  then  used  in  equation  (70)  to  calculate  a  in  order  to  detennine  the  equivalent  thickness 
and  stiffness  E*  and  t*. 


b* 


Eskin 


(83) 


Where: 

b  —  stiffener  thickness 
b*  —  transformed  stiffener  thickness 
Estif  —  stiffener  modulus  of  elasticity 
Eskin  —  skin  modulus  of  elasticity 

After  transfonning  the  stiffener  thickness  it  is  possible  to  treat  the  entire  structure 
as  having  isotropic  properties  of  constant  stiffness  in  order  to  use  equation  (66).  This 
process  was  used  to  determine  the  W/B  of  an  isogrid  stiffened  sphere  assuming  UHM 
carbon  as  described  above.  The  resulting  W/B  is  displayed  in  Figure  38.  The  W/B  is 
slightly  better  than  beryllium,  but  LTA  flight  is  not  achieved.  Note  that  the  bottom  of  the 
curve  in  Figure  38  is  curved,  not  flat  like  in  Figure  37.  This  is  because  the  compressive 
strength  of  the  epoxy  matrix,  as  shown  in  Table  7,  is  sufficient  to  prevent  material  failure 
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and  allow  very  small  t.  In  this  case  the  E/p2  value  is  the  active  material  constraint,  and 
changes  in  stiffness  would  result  in  smaller  W/B. 


Figure  38:  UHM  Carbon  Epoxy  Isogrid  Optimization 


The  optimal  design  values  for  a  sphere  of  rs  =  0.33  m  are  as  follows:  h  — 

2.1  mm,  t  —  0.042  mm.,  b  —  0.035  mm,  and  d  —  1  mm.  This  design  would  be 
extremely  difficult  to  build  on  account  of  the  small  size  and  high  number  of  stiffeners 
required.  In  order  to  arrange  triangular  stiffeners  on  the  surface  of  a  sphere,  a  geodesic 
arrangement  must  be  used.  Therefore  the  geometry  used  to  calculate  the  geodesic  sphere 
can  also  be  used  to  calculate  the  arrangement  of  isogrid  stiffeners.  An  rs  —  0.33  m  with 
an  a  =  2.4  mm  implies  that  the  major  triangle  of  length  0.347  (where  a  =  1.05  ■  rs), 
would  be  divided  144  times  to  get  a  minor  triangle  length  of  4.1  mm,  meaning  /=144. 
Using  equation  (28),  the  number  of  individual  blade  type  stiffeners  (triangle  legs)  would 
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be  622,080.  The  Isogrid  Handbook  [36]  actually  mentions  that  a  feasible  way  to  apply  an 
isogrid  to  a  sphere  is  using  the  geodesic  method  described  in  [37].  A  close-up  view  of  the 
surface  of  the  proposed  sphere  with  f=85  is  shown  in  Figure  39.  An  f=  1 44  would  be  even 
more  complex.  A  feasible  method  for  building  such  a  complex  structure  would  require 
automation. 


Figure  39:  Isogrid  Visualization  /=85 


After  conducting  the  optimization  technique  of  simultaneously  satisfying  all  three 

constraints  it  became  apparent  that  the  isogrid  did  not  perfonn  as  well  as  expected.  The 

possible  cause  was  that  the  optimal  design  did  not  have  all  constraints  active.  Finding  this 

point  required  a  more  robust  optimization  method.  Therefore  the  isogrid  was  optimized 

using  the  non-linear  programming  optimization  routine  called  fmincon.  The  program  was 

inconsistent  with  its  results  at  first  because  of  scaling  problems.  The  variables  h,  tskiTls , 

b,  and  d  were  so  small  that  the  program’s  step  size  tended  to  overshoot  the  optimal  point 

and  failed  to  converge  on  the  right  value.  After  properly  scaling  the  problem  to  make  the 

variables  have  a  magnitude  close  to  unity,  fmincon  consistently  converged  on  a 

106 


W/B=0.81.  This  occurred  at  h=\A  mm,  tskin  = 0.032  mm,  6=0.015  mm,  and  d=\21  mm. 
These  values  did  not  lie  in  the  two  dimensional  design  space  shown  in  Figure  38  because 
not  all  constraints  were  active  at  the  optimal  point.  Only  the  global  buckling  constraint, 
equation  (66),  was  active  with  a  value  of  -0.01.  The  negative  value  indicates  that  the 
constraint  was  satisfied,  the  magnitude  indicates  it  was  just  barely  satisfied,  thus  it  was  an 
active  constraint.  The  other  constraints  were  inactive  with  values  for  the  skin  buckling  of 
-63  and  rib  buckling  of  -73.  The  negative  values  indicate  the  constraints  were  satisfied, 
which  means  the  buckling  load  was  not  exceeded  by  the  stress.  The  magnitudes  indicate 
that  the  skin  and  stiffeners  had  thicker  cross  sections  than  the  minimum  required  to 
prevent  buckling  of  the  skin  or  stiffeners,  which  means  these  constraints  were  not  active. 

Two  additional  constraints  were  added  to  the  optimization.  These  were  the 
requirement  that  the  skin  compressive  yield  strength  and  stiffener  compressive  yield 
strength  not  be  exceeded  by  the  stress  in  the  structure.  This  resulted  in  a  total  of  5 
equations  and  four  variables.  The  problem  was  not  over-constrained,  however,  as  long  as 
no  more  than  four  variables  were  active.  It  turned  out  that  at  the  optimal  point  for  each  of 
the  isogrid  designs  only  one  constraint  was  ever  active,  the  global  buckling  constraint, 
equation  (66). 

The  same  optimization  was  run  for  pure  beryllium  resulting  in  a  W/B=l  .7,  just  as 
in  the  method  of  simultaneously  solving  constraints.  This  was  the  expected  result  because 
the  compressive  strength  of  beryllium  prevented  the  design  from  falling  below  a 
minimum  cross  sectional  area,  which  prevented  reduction  in  W/B  using  any  method  of 
optimization.  No  particular  result  was  found  because  there  were  a  range  of  values  for 
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which  a  W/B=1.7  could  occur,  just  as  shown  in  Figure  37  using  the  first  optimization 
method. 

The  same  optimization  was  run  for  a  sphere  with  a  beryllium  skin  and  a  carbon 
fiber  stiffener.  A  W/B=0.79  was  obtained  with  a  geometry  of  h= 1.3  mm,  tskiris=0.026 
mm,  6=0.014  mm,  and  <7=1.2  mm.  This  hybrid  structure  of  beryllium  and  UHM  carbon 
composite  combined  the  higher  isotropic  stiffness  of  Beryllium  with  the  higher  axial 
stiffness,  lower  density,  and  higher  compressive  strength  of  UHM  carbon. 

The  geometry  and  W/B  of  each  isogrid  design  is  listed  in  Table  9.  Note  the 
comparatively  large  value  for  b  and  small  value  for  d  for  the  beryllium  skin/I  JHM  carbon 
stiffener  hybrid  structure.  This  can  be  explained  by  the  high  stiffness  in  the  beryllium 
skin  which  did  not  require  a  substantial  stiffener  to  be  added  to  prevent  skin  or  global 
buckling.  The  skin  did  require  additional  area,  however,  to  spread  out  the  stress  due  to  its 
low  compressive  yield  stress.  The  UHM  carbon  epoxy  is  substantially  less  dense  than  the 
beryllium,  which  makes  it  a  better  material  for  adding  thickness. 
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Table  9:  Isogrid  Results 


Material 

Optimization 

Method 

h 

t 

b 

d 

W/B 

Diagram 

Beryllium 

Simultaneous 

Constraint 

Satisfaction  Non¬ 
linear 

Programming 

1.8E- 

03 

3.3E- 

05 

4.8E- 

05 

1.1E- 

03 

1.66 

UHM 

Carbon 

Epoxy 

Simultaneous 

Constraint 

Satisfaction 

2.1E- 

03 

4.2E- 

05 

3.5E- 

05 

1.0E- 

03 

1.03 

UHM 

Carbon 

Epoxy 

Non-linear 

Programming 

1.4E- 

03 

3.2E- 

05 

1.5E- 

05 

1.3E- 

03 

0.81 

Beryllium 
Skin  UHM 

Carbon 

Epoxy 

Stiffeners 

Non-linear 

Programming 

1.3E- 

03 

2.6E- 

05 

4.8E- 

05 

1.4E- 

05 

0.79 

- 

4.5  Geodesic  Sphere 

The  geodesic  sphere  results  are  given  in  this  section,  as  well  as  a  description  of 
the  design  space.  The  geodesic  sphere  was  evaluated  for/=l  through  10.  If  the 
optimization  showed  that  higher / resulted  in  lower  W/B,  then  higher  values  of f  would 
have  been  investigated.  This  approach  was  taken  because  of  the  high  computational  cost 
of  generating  and  evaluating  structures  with  high  numbers  of  elements.  The  number  of 
elements  in  a  geodesic  structure  was  equal  to  the  number  of  edges,  ne,  given  by  equation 
(30).  The  computational  cost  of  the  optimization  routine  depended  on  f  with /=  1  (30 
elements)  having  an  average  run  time  of  1  second  and  an /= 9  (2430  elements)  having  an 
average  run  time  of  120  seconds  per  iteration  on  an  Intel®  Core  ™  i7-2630QM  CPU  @ 
2.00  GHz  with  8  MB  of  RAM.  This  trend  is  plotted  in  Figure  40.  The  differences  in  run 
time  were  primarily  due  to  the  FEA  portion  of  the  routine  with  the  large  number  of 


109 


elements  at  higher  frequencies  requiring  significant  computation  times.  The  graph  in 
Figure  40  should  be  quadratic,  which  it  appears  to  be,  Because  ne  is  a  function  of  f2. 


Figure  40:  Effect  of  Geometric  Frequency  on  Run  Time 


It  turned  out  that  lower  frequencies  resulted  in  lower  W/B,  so  that  no/>10  was 
investigated.  There  are  two  reasons  an/=  1  is  optimal.  The  first  is  minimal  structure.  At 
higher  frequencies  the  sphere  approaches  a  solid  sphere,  whereas  low  frequencies  have 
more  empty  space.  This  can  be  seen  in  Figure  41,  where  the  isocost  lines  (lines  of 
constant  cost)  which  correspond  to  geodesic  frame  W/B,  move  left  with  increasing/  The 
axes  of  Figure  41  are  the  variables  rCb  and  tcb-  For  a  given  combination  of  rCb  and  tCb,  an 
increase  in/ results  in  an  increase  in  W/B.  This  makes  sense  because  it  is  equivalent  to 
keeping  the  cross  section  of  the  cylindrical  beams  constant  and  increasing  total  length  of 
the  beams  (if  all  members  were  added  together  to  fonn  one  long  cylinder).  This  would 
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increase  overall  weight  of  the  structure  while  keeping  buoyancy  constant,  which  would 
increase  W/B. 


Figure  41:  Effect  of  Geometric  Frequency  on  Feasible  Design  Space 


The  other  reason /=  1  is  optimal  is  that  the  icosahedron  is  perfectly  symmetric,  so 

that  every  beam  has  the  same  stress  when  loaded  evenly.  Figure  42  shows  the 

relationship  between  stress  and  f  At  higher/ the  average  stress  in  all  the  beams  decreases 

in  magnitude,  however  some  members  have  more  stress  than  others  due  to  asymmetrical 

loading.  The  icosahedron  (/=1)  is  the  structure  with  the  lowest  maximum  stress  in  its 

members  because  it  is  a  perfectly  symmetric  structure,  so  that  its  maximum  stress 

magnitude  is  equal  to  its  average  stress  magnitude.  This  property  of  the  icosahedron,  in 
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addition  to  having  the  lowest  structural  mass,  results  in  the  lowest  W/B  of  any  geodesic 
structure  investigated.  Note  that  Figure  42  is  a  plot  of  compressive  stress,  so  that  all 
values  are  negative  and  a  higher  compressive  stress  magnitude  is  lower  on  the  y  axis  than 
a  lower  compressive  stress  magnitude.  The  solid  dots  indicate  maximum  stress  anywhere 
in  the  structure,  while  the  hollow  dots  indicate  mean  stress  throughout  the  structure.  For 
the  icosahedron,  which  corresponds  to  a  geodesic  sphere  with  an/=l,  the  mean  stress 
equals  the  maximum  stress. 


Figure  42:  Effect  of  Geometric  Frequency  on  a 
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Figure  43  shows  the  design  space  for/=l.  White  denotes  the  feasible  design  space 
with  isocost  lines  denoting  W/Bs  of  1  through  4.  The  constraints  gi,  g2  and  g3  are  labeled 
with  the  infeasible  design  space  shown  in  color.  The  constraint  g4  is  not  pictured  due  to 


112 


being  off  the  right  side  of  the  graph.  The  optimal  point  shows  a  W/B  less  than  1 .  The  gi 
constraint  is  pictured  in  red  (brown  when  it  overlaps  go).  The  gi  constraint  is  the 
requirement  for  compressive  yield  strength  of  the  material  to  exceed  the  maximum 
compressive  stress  in  the  structure.  As  would  be  expected,  this  constraint  is  parallel  to  the 
isocost  lines  because  cross  sectional  area  is  directly  related  to  both  weight  and  stress.  If 
cross  sectional  area  is  always  exactly  sufficient  to  satisfy  gi,  it  will  correspond  to  the 
same  amount  of  material  in  the  frame,  which  will  correspond  to  a  constant  weight.  Cross 
sectional  area  is  a  product  of  the  two  variables  rcb  and  tCb,  which  are  the  x  and  y  axes  of 
Figure  43.  The  advantage  of  having  gi  be  the  active  constraint  is  that  it  allows  for 
flexibility  in  design.  An  rcb  anywhere  between  4mm  and  1cm,  as  shown  on  Figure  43, 
will  result  in  a  W/B  for  the  frame  of  0.57.  If,  for  example,  the  minimum  thickness  that 
can  be  constructed  of  a  pultruded  UHM  carbon  epoxy  rod  happens  to  be  0. 15  mm,  then 
that  thickness  can  be  selected,  and  the  corresponding  rcb=5  mm  could  be  used  without  an 
increase  in  W/B. 

The  requirement  for  the  buckling  load  of  the  cylindrical  beams  to  exceed  the  load 
on  any  member  in  the  structure  is  defined  by  go.  The  constraint  go  is  optimal  at  infinite 
radius  and  infinitesimal  thickness.  As  can  be  seen  in  Figure  43,  minimum  rcb  is  defined 
by  go,  and  minimum  tcb  is  defined  by  gi. 

The  constraint  g3  is  the  requirement  that  thickness  not  exceed  half  the  radius,  and 
so  is  only  active  for  high  thickness  and  low  radius.  This  can  be  seen  on  Figure  43  where 
g3  is  only  active  at  low  radii  and  high  thickness. 

The  constraint  g4  does  not  even  appear  on  Figure  43.  This  is  because  g4  sets  an 

upper  limit  on  rcb  as  a  function  of f.  Higher  fs  result  in  more  restrictions  on  rcb,  so  that  g4 
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only  becomes  a  factor  in  design  spaces  with  high /  The  constraint  g4  is  visible  in  Figure 
41  on  the  graph  of f=5. 
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Figure  43:  Graphical  Representation  of  Geodesic  Sphere  Design  Space 


Radius,  meters 


A  minimum  W/B=  0.57  was  attained  for  the  UHM  carbon  epoxy  frame  with/=l, 
rch  —  1  cm,  and  tcb  =  0.64  mm  for  the  frame  alone  as  shown  in  Figure  44.  The  cross 
section  of  a  cylindrical  tube  is  shown  with  the  thickness  of  the  cylinder  wall  scaled 
proportional  to  the  radius  of  the  cylinder.  The  icosahedron  is  also  shown  with  the  line 
thickness  scaled  proportional  to  the  triangle  leg  length.  The  fmincon  program  proved 
robust  enough  for  the  problem  even  though  it  was  not  designed  to  handle  discrete  design 
variables.  Although / was  a  discrete  design  variable,  fmincon  was  usable  by  rounding/ to 
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the  nearest  integer  within  the  subroutines  of  the  optimization  routine,  as  well  as  rounding 
the  final  value  of f  at  the  end  of  the  optimization  prior  to  evaluating  it  in  the  cost  function. 


Figure  44:  Optimal  Geometry  for  rs=0.33  m  Geodesic  Sphere 


The  skin  thickness,  tsk,  was  not  optimized  but  was  assumed  to  be  5  mm  thick 
Vectran  for  the  purpose  of  FEA  evaluation  of  a  triangular  plate  under  1.5  times  sea  level 
pressure  with  clamped  edges  as  shown  in  Figure  45  .  It  was  found  that  with  a  thickness  of 
5mm,  a  skin  with  Vectran  material  properties  as  shown  in  Table  6,  that  the  maximum 
displacement  at  the  center  of  the  skin  was  0.59  mm.  The  maximum  stress  in  the  skin  was 
49  MPa  which  is  far  below  the  maximum  strength  of  Vectran  of  3.2  GPa  as  shown  in 
Table  6.  This  only  occurred  on  the  edges  and  center  of  the  triangle,  whereas  the  average 
stress  in  the  skin  was  below  30  MPa.  An  optimized  skin  could  be  reinforced  in  the  places 
of  greatest  stress,  so  that  the  average  thickness  could  be  designed  for  a  stress  of  30  MPa. 
One  of  the  highest  specific  tensile  strengths  in  available  materials  is  found  in  Zylon,  with 

a  tensile  strength  of  5.8  GPa  [33].  Assuming  the  skin  thickness  of  5  mm  and  a  stress  of  30 
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MPa  from  the  FEA  shown  in  Figure  45,  the  force/unit  length  was  calculated  in  the 


structure  by  multiplying  tskin  by  stress.  Taking  this  force/length  and  dividing  by  5.8  GPa 
gave  a  minimum  skin  thickness  of  0.0259  mm.  This  process  is  shown  in  equation  (84). 


aFEA  ‘ 


tskinFEA 


(J-\, 


=  t 


skin 


required. 


(84) 


Where: 

tskin.pE A  —  skin  thickness  used  in  FEA 

tskin  ^  =  minimum  required  skin  thickness 

required  1 

(jFEa  =  Stress  calculated  by  FEA 
oy  =  yield  strength  of  material 

The  minimum  required  skin  thickness,  tskinrequired  =  0.0259  mm,  gave  a 
W/B=0.37  for  the  skin  and  a  W/B=0.94  for  the  whole  structure.  The  payload  capability  of 
a  structure  with  a  W/B=0.94  is  shown  in  Figure  46.  The  mass  of  the  structure  would  be 
0.94/(l-.94)=15.7  times  the  payload.  For  a  diameter  of  1  meter,  the  structure  would  have 
a  mass  of  580  grams  and  a  payload  of  37  grams. 
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Figure  45:  FEA  of  Vectran  Triangular  Plate  With  Clamped  Edges 
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Figure  46:  Payload  Capability  for  W/B=0.94 


A  much  more  promising  material  to  be  used  for  the  membrane  is  graphene. 
Graphene  has  a  tensile  strength  of  135  GPa,  [58],  which  would  enable  a  skin  only  0.001 
mm  thick.  A  skin  of  graphene  would  have  a  W/B=0.001,  which  would  result  in  an  overall 
W/B=0.57.  This  would  be  an  effectively  weightless  skin.  Currently,  graphene  appears  to 
only  be  available  in  small  sizes  on  the  order  of  several  centimeters  in  length,  but 
graphene  represents  the  future  potential  of  material  research  [58]. 

Optimization  of  the  skin  would  be  a  tradeoff  between  loss  of  buoyancy  due  to 
decreasing  volume  (  Vgs)  with  skin  displacement  as  a  function  of  tsk,  and  a  decrease  in 
weight  as  a  function  of  tsk .  The  loss  in  Vgs  would  occur  when  the  skin  displaced  inward 
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due  to  external  pressure,  thus  reducing  the  total  volume  of  displaced  air  contributing  to 
buoyancy.  Diffusivity  is  a  measure  of  the  rate  at  which  air  could  seep  through  the  Mylar 
membrane.  If  this  occurred  the  weight  of  the  structure  would  increase  due  to  added  air 
molecules  inside  the  skin  of  the  structure.  This  is  the  same  process  as  a  leaking  ship  that 
sinks  as  it  takes  on  water.  Optimization  of  both  W/B  and  diffusivity  would  be  a  multi¬ 
objective  problem  in  which  the  diffusivity  of  air  through  the  reinforced  Mylar  membrane 
as  a  function  of  tsk  would  be  considered.  The  cost  function  would  need  to  be  adjusted  to 
incorporate  diffusivity  with  weighting  values  assigned  to  both  W/B  and  diffusivity.  A 
Pareto  front  could  then  be  constructed  to  give  designers  a  choice  between  design  goals.  A 
weighting  function  could  also  be  fonned  such  as  maximizing  mission  time,  where 
mission  time  is  a  function  of  both  diffusivity  and  W/B.  A  Pareto  front  is  constructed  by 
giving  different  weights  to  the  objective  functions  of  a  multi-objective  optimization 
problem.  If  one  objective  is  given  a  lot  more  weight  than  the  other,  the  problem  will  be 
optimized  primarily  for  that  objective  alone.  If  both  objectives  are  given  equal  weight  a 
solution  that  is  half  way  between  optimizing  either  objective  will  be  the  optimal  point. 

4.6  Potential  Missions 

The  vacuum  LTA  vehicle  has  the  potential  to  fulfill  missions  in  transportation, 
surveillance,  communications  relay,  and  sensor  platforms  just  like  any  other  LTA 
vehicle.  Whether  a  vacuum  LTA  vehicle  is  suited  for  any  of  these  missions  will  depend 
on  each  of  the  design  challenges  covered  previously.  The  most  important  hurdle  the 
vacuum  structure  must  clear  is  becoming  LTA.  Therefore  minimizing  W/B  is  the  first 
priority.  Based  on  W/B,  a  von  Karman  Efficiency  can  be  predicted  for  a  theoretical 
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vehicle  that  uses  a  vacuum  LTA  structure.  This  von  Kannan  efficiency  in  turn  predicts 
whether  the  vehicle  could  be  economically  used  for  transportation. 

One  of  the  questions  posed  early  in  this  study  was  if  there  was  an  ideal  altitude  at 
which  a  vacuum  LTA  structure  could  be  designed  to  operate.  Figure  47  is  a  plot  of  W/B 
for  a  vacuum  LTA  structure  at  various  altitudes  using  altitude  tables  for  a  standard  day 
[28].  Figure  33  and  Figure  47  were  created  using  the  Matlab  program  titled  LTApvalt.m  , 
Appendix  A.  1 1 .  There  are  two  lines  with  positive  slope  in  Figure  47.  The  solid  line 
represents  a  vacuum  LTA  structure  designed  to  achieve  W/B=0.95  at  sea  level.  At  higher 
altitudes,  buoyancy  decreases  due  to  decreased  air  density,  and  therefore  W/B  increases. 
The  dashed  line  represents  the  minimum  W/B  for  vacuum  LTA  structures  specifically 
designed  to  resist  the  pressure  of  a  particular  altitude  (and  no  more).  The  dash  dot  line 
represents  a  vacuum  LTA  structure  designed  for  500  meters  pressure  altitude.  This 
structure  would  have  a  lower  W/B  than  the  vacuum  LTA  structure  designed  for  sea  level, 
but  would  have  a  minimum  altitude  of  500  meters.  Descending  below  500  meters 
pressure  altitude  would  cause  the  structure  to  exceed  its  design  pressures,  which  could 
cause  it  to  fail.  Even  though  the  structure  designed  for  500  meters  has  a  better  W/B  at 
higher  altitudes  than  the  structure  designed  for  sea  level,  the  minimum  W/B  achievable  is 
by  the  structure  designed  for  sea  level,  while  it  is  at  sea  level.  This  study  shows  that  the 
optimal  elevation  for  W/B  is  always  the  lowest  available  elevation.  Therefore  the  lowest 
W/B  achievable  by  a  vacuum  LTA  structure  is  at  sea  level.  For  this  reason  all  structures 
evaluated  during  the  course  of  this  study  were  evaluated  for  the  W/B  achievable  at  sea 
level  standard  atmosphere  conditions  in  order  to  give  a  best  case  scenario  as  well  as  for 
the  purpose  of  continuity. 
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Figure  47:  Effect  of  Increasing  Designed  Deployment  Elevation 


Figure  48  shows  the  W/B  for  LTA  structures  of  various  W/Bs  at  sea  level.  In 
order  for  a  LTA  structure  to  operate  at  high  altitudes,  it  must  have  a  low  W/B  at  sea  level. 
The  maximum  altitude  for  a  given  LTA  structure  is  the  altitude  at  which  its  W/B=l, 
which  is  the  shaded  region  on  both  Figure  47  and  Figure  48.  This  is  the  main  reason 
cargo  LTA  vehicles  were  commonly  operated  below  1000  feet  [9].  Sea  level  provides  the 
maximum  possible  buoyancy  for  any  LTA  vehicle  and  is  therefore  the  ideal  place  to 
operate  cargo  oriented  LTA  vehicles.  Without  low  W/B,  any  mission  requiring  high 
altitude  will  be  impossible  for  the  vacuum  LTA  vehicle.  A  W/B=0.09  would  be  required 
for  a  LTA  vehicle  to  operate  at  the  design  altitude  of  the  Lockheed  HALED  [59]. 
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Figure  48:  Effect  of  Altitude  on  W/B 
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The  von  Karman  efficiency  of  the  design  space  for  a  vacuum  LTA  vehicle  is 
graphed  in  Figure  50.  This  vehicle  was  assumed  to  have  a  W/B=  0.8,  an  aerodynamic 
shroud,  and  a  BLC  device  added  that  gave  the  same  CD  —  0.0169  as  the  BLC  airship 
previously  evaluated.  The  assumption  was  made  that  the  shroud  would  increase  the  W/B 
to  0.9  and  that  it  would  give  the  spherical  vacuum  structure  a  more  aerodynamic  shape.  A 
conceptual  design  of  this  vehicle  is  shown  in  Figure  49  and  is  composed  of  the  Griffith 
Airfoil  used  as  a  body  of  revolution  as  described  and  drawn  by  Goldschmied  [24],  with 
three  vacuum  LTA  spheres  inside. 
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Figure  49:  Conceptual  Design  of  BLC  Vacuum  LTA  Vehicle  [24] 
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Von  Karman  efficiency  was  plotted  for  sea  level  standard  day  conditions  and 
payload  was  calculated  according  to  equation  (4).  The  lift  for  this  vehicle  was  abysmal 
but  the  von  Karman  efficiency  was  competitive  due  to  low  drag.  If  the  drag  were 
increased  to  the  value  of  a  smooth  sphere  where  CD  —  0.1  [5 1]  and  the  W/B  remained  at 
0.8,  the  resulting  von  Karman  efficiency  of  the  design  space  would  look  like  Figure  5 1 . 
This  shows  that  the  increased  drag  has  a  significant  effect  on  von  Karman  efficiency, 
more  so  than  the  ten  percent  change  in  W/B.  Therefore  as  long  as  the  shroud  could  be 
made  light  enough  for  the  structure  to  float,  there  is  a  good  chance  that  improving  the 
aerodynamics  of  the  vacuum  LTA  structure  beyond  that  of  a  sphere  would  be  beneficial. 
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Von  Karman  Efficiency 


Figure  50:  von  Karman  Efficiency  at  W/B=0.8  CD  —  0.  0169 
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Figure  51:  von  Karman  Efficiency  at  W/B=0.8  CD  =  0. 1 


The  von  Karman  efficiency  of  a  vacuum  LTA  vehicle  with  a  W/B  =0.8  or  higher 
at  sea  level  is  unlikely  to  be  used  for  transporting  cargo.  This  conclusion  was  made  based 
on  the  altitude  limitations  shown  in  Figure  48  and  the  poor  von  Karman  efficiency  shown 
in  Figure  50  and  Figure  5 1 .  If  the  BLC  can  be  achieved  on  the  vacuum  LTA  vehicle  it 
would  perfonn  a  little  worse  than  the  Zeppelin  NT  for  a  volume  equal  to  the  Zeppelin 
NT.  Any  BLC  that  could  be  applied  to  the  vacuum  LTA  vehicle  could  also  be  applied  to 
a  LTA  vehicle  relying  on  a  lifting  gas,  so  that  the  vacuum  LTA  vehicle  is  at  an  efficiency 
disadvantage  to  the  lifting  gas  airship.  The  vacuum  LTA  vehicle  does  pose  advantages, 
however,  due  to  its  lack  of  need  for  a  lifting  gas,  which  would  enable  its  use  in  areas 
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where  shipping  helium  containers  is  cost  prohibitive,  or  where  the  danger  of  fire  due  to 
hydrogen  is  unacceptable.  Therefore  it  is  conceivable  that  a  mission  suited  for  the 
vacuum  LTA  vehicle  would  be  one  requiring  low  payload  and  low  altitude,  such  as  an 
Intelligence,  Surveillance,  Reconnaissance  (ISR)  sensor  platform  for  use  in  the  urban 
environment.  These  missions  could  be  conducted  in  hostile  environments  such  as  around 
chemical,  biological,  or  radiation  hazards.  A  geodesic  sphere  could  be  constructed  so  that 
certain  beams  were  designed  to  buckle  without  damage.  This  would  allow  intentional 
buckling  of  the  structure  for  storage  purposes.  When  the  structure  needed  to  be  re¬ 
deployed  it  could  be  inflated  to  unbuckle  the  beams,  deflated  to  cause  positive  buoyancy, 
and  re-deployed.  This  would  require  a  valve  and  a  vacuum  pump,  which  would  require  a 
power  source,  but  not  a  lifting  gas. 

A  vacuum  LTA  vehicle  with  a  radius=0.33  meters  could  fit  through  a  standard 
door.  If  it  had  a  W/B=0.8,  it  would  have  a  total  structure  mass  of  180  grams  and  be  able 
to  carry  46  grams  at  sea  level.  Any  control  systems,  fuel,  motors  and  payload  would  have 
to  weigh  less  than  46  grams. 

4.7  Summary 

The  rotating  cylinder  vehicle  was  deemed  infeasible  due  to  its  high  energy 
requirements  due  to  torque  drag.  Even  a  rotating  cylinder  vehicle  with  a  drag  torque 
reduction  of  99%  below  a  single  rotating  cylinder  would  be  infeasible  at  sea  level.  At  20 
km  MSL  the  RCVLTAV  with  a  99%  drag  torque  reduction  would  appear  energetically 
viable,  however  the  requirement  to  be  lifted  to  altitude  by  an  external  aircraft,  and  the 
requirement  for  a  99%  reduced  drag  torque  make  the  RCVLTAV  appear  infeasible. 
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A  grid  stiffened  sphere  was  optimized  using  two  methods,  one  in  which  all 
constraints  were  simultaneously  satisfied,  and  a  non-linear  programming  optimization 
problem.  Two  material  arrangements  were  found  with  a  W/B<1.  The  first  was  an  UHM 
carbon  epoxy  sphere  with  unidirectional  stiffeners  and  a  quasi-isotropic  skin  which  was 
predicted  to  achieve  a  W/B=0.81.  The  second  was  a  beryllium  skin  with  UHM  carbon 
fiber  unidirectional  stiffeners  which  achieved  a  W/B=0.79.  This  was  achieved  with  a  high 
geometric  frequency  of  very  thin  stiffeners,  which  calls  into  question  the  feasibility  of 
construction.  Advanced  automated  methods  of  composite  construction  would  be  required 
to  construct  such  high  frequency  geodesic  grids. 

A  geodesic  frame  composed  of  cylindrical  beams  made  of  UHM  carbon  fiber  was 
predicted  to  achieve  a  W/B=0.57.  A  skin  made  of  Mylar  reinforced  with  Zylon  was 
predicted  to  add  an  additional  W/B=0.37,  which  would  result  in  an  overall  W/B=0.94.  A 
skin  with  the  properties  of  graphene  would  have  a  W/B=0.001,  which  would  result  in  a 
vehicle  with  a  W/B=0.57.  This  conclusion  was  attained  using  FEA  in  conjunction  with 
Euler  buckling  loads  and  published  material  properties.  Optimization  techniques  included 
graphical,  the  genetic  algorithm  function  in  Matlab,  and  the  fmincon  function  in  Matlab. 

The  von  Karman  efficiency  for  a  vacuum  LTA  vehicle  with  a  structural  W/B=0.8 
was  evaluated.  The  vacuum  LTA  vehicle  could  not  compete  with  a  lifting  gas  LTA 
vehicle  as  far  as  efficiency  is  concerned.  It  would,  however,  have  the  unique  property  of 
not  relying  on  a  lifting  gas.  This  could  make  the  vacuum  LTA  vehicle  suitable  for 
missions  in  which  efficiency  was  not  of  primary  concern. 
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V.  Conclusions  and  Recommendations 


5.1  Chapter  Overview 

This  study  began  by  investigating  LTA  vehicles  to  include  their  structure, 
missions,  and  design  challenges.  The  buckling  of  vacuum  structures  was  combined  with 
buoyancy  to  find  several  structures  that  could  possibly  achieve  positive  buoyancy  using  a 
vacuum.  The  performance  of  vehicles  using  these  vacuum  structures  were  then  predicted 
using  von  Karman  Efficiency. 


5.2  Conclusions  of  Research 

Eight  questions  were  posed  in  the  introduction  section.  The  answers  to  these 
questions  are  available  in  Table  10. 


Table  10:  Research  Questions  Revisited 


1 

Can  positive  buoyancy  be  achieved  by  a  structure  in  air  using  a  vacuum? 

Probably.  The  grid  stiffened  structure  could  possibly  achieve  positive 
buoyancy  using  a  vacuum  but  advanced  automated  construction  of  large 
numbers  of  stiffeners  would  be  required.  The  beryllium  and  UHM 
carbon  epoxy  hybrid  was  the  most  successful  structure  investigated,  with 
a  predicted  W/B=0.79.  The  geodesic  sphere  appears  to  provide  a  feasible 
means  by  which  positive  buoyancy  can  be  achieved  by  a  vacuum 
structure  using  available  materials  such  as  UHM  carbon  tubes  and  Zylon 
reinforced  skin.  A  W/B=0.94  for  the  geodesic  structure  was  predicted, 
with  W/B=0.57  attributed  to  the  frame,  and  W/B=0.37  attributed  to  the 
skin.  Graphene  would  far  exceed  the  material  requirements.  The  rotating 
cylinder  vehicle  cannot  feasibly  achieve  positive  buoyancy. 

2 

Is  a  vacuum  LTA  structure  a  feasible  platform  for  any  existing  missions? 

Yes.  The  vacuum  LTA  structure  could  potentially  be  used  for  any 
mission  not  requiring  a  low  W/B.  This  would  include  any  mission 
suitable  for  the  characteristics  of  a  LTA  vehicle  that  do  not  require  a  low 
von  Kannan  efficiency.  A  mission  that  appears  particularly  suited  for  the 
vacuum  LTA  vehicle  is  a  low  altitude  sensor  platform  for  ISR  in  hostile 
environments. 

3 

What  is  the  ideal  internal  pressure  of  a  partial  vacuum  LTA  structure? 

The  closer  the  internal  pressure  can  be  to  absolute  vacuum,  the  better 

128 


W/B  achievable  by  the  vacuum  LTA  structure. 

4 

What  is  the  ideal  operating  altitude  of  a  vacuum  LTA  structure? 

Sea  level  is  the  ideal  operating  altitude.  Any  increase  in  altitude  results 
in  an  increase  in  W/B  as  shown  in  Figure  47. 

5 

Are  there  viable  geometric  shapes  for  a  vacuum  LTA  structure? 

The  sphere  is  the  best  geometric  shape,  but  geodesic  approximations  to 
the  sphere  appear  viable  as  well.  The  icosahedron  actually  achieved 
lower  W/B  than  more  complicated  geodesic  shapes  due  to  its  symmetry 
and  sparse  structure. 

6 

Do  construction  techniques  exist  which  can  construct  the  required 
geometric  shapes? 

The  answer  appears  unlikely  for  the  grid  stiffened  structure  as  long  as  it 
is  made  of  carbon  fibers.  It  is  possible  that  a  dissolvable  mandrill  could 
be  constructed  by  a  rapid  prototype  machine  and  carbon  fiber  applied 
using  continuously  wound  fibers  using  an  automated  filament  winder. 

The  thicknesses  of  the  skin  and  stiffeners  of  the  grid  stiffened  structure 
would  also  pose  a  problem  unless  the  structure  was  made  very  large.  It 
appears  the  geodesic  sphere  could  be  constructed  because  the  optimal 
shape  is  composed  of  only  30  tubes  of  reasonable  thickness.  The 
challenge  to  building  this  structure  lies  in  connecting  the  tubes  at  the 
vertices. 

7 

What  are  the  critical  material  properties  for  a  vacuum  LTA  structure? 

E  —  Young's  modulus  of  elasticity 
acy  —  compressive  yield  strength 
p  =  material  density 
oy  =  tensile  yield  strength 

A  material  search  can  be  perfonned  for  the  geodesic  frame  or  isogrid 
structure  by  looking  for  the  maximum  E / p2 .  The  ocy  must  also  be 
checked  as  it  may  be  the  active  constraint  as  in  the  case  of  Beryllium. 

For  skin,  permeability  is  important,  as  well  as  tensile  yield  strength  oy 
and  p.  A  material  search  for  the  ideal  skin  could  be  found  by  searching 
for  the  maximum  <Jy/p 

8 

Do  materials  exist  from  which  a  vacuum  LTA  structure  can  be 
constructed? 

Yes.  UHM  carbon  fiber  in  an  epoxy  matrix  is  the  most  likely  material  to 
be  used  for  the  structure.  Stiffer  fibers  and  less  dense  matrix  materials 
are  the  most  promising  area  for  material  improvement.  For  the 
membrane  a  high  specific  strength  material  like  Zylon  appears  most 
promising. 
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5.3  Significance  of  Research 

A  vacuum  LTA  structure  would  replace  the  need  for  a  lifting  gas.  If  the  risk  of 
fire  is  unacceptable  for  a  particular  mission  and  helium  is  unavailable  due  to  reduced 
supply,  vacuum  LTA  structures  may  be  the  best  alternative.  At  the  current  state  of 
research  a  vacuum  structure  can  not  achieve  W/Bs  as  low  as  LTA  structures  containing 
helium  or  hydrogen  lifting  gases.  Nevertheless,  the  construction  of  a  vacuum  LTA 
structure  would  constitute  a  technological  milestone  in  structural  design.  The  largest 
contribution  of  this  research  is  that  it  appears  to  be  the  first  time  FEA  and  non-linear 
programming  optimization  techniques  have  been  used  to  show  that  a  structure  can 
achieve  positive  buoyancy  in  air  using  a  vacuum. 

5.4  Recommendations  for  Future  Research 

Recommended  future  research  includes  optimization  of  the  skin  for  diffusivity 
and  deflection,  evaluation  of  the  bending  effects  due  to  membrane  loads  on  the  sides  of 
the  beams  using  FEA,  and  experimental  analysis  of  both  membrane  deflection  and  global 
buckling  of  the  icosahedron  in  order  to  validate  the  FEA  models.  An  experiment  could  be 
perfonned  on  inexpensive  materials  such  as  fiberglass  in  order  to  perfect  construction 
techniques  and  validate  W/B  predictions.  This  experiment  could  be  performed  by 
building  the  structure  with  geometry  optimized  for  the  cheaper  material,  predicting  the 
pressure  at  which  it  will  fail  using  the  cheaper  material,  then  removing  the  air  from  inside 
of  the  structure  using  a  vacuum  pump  while  measuring  the  pressure  inside  and  outside  of 
the  structure  to  observed  the  critical  pressure  at  which  the  structure  fails.  This  would  be  a 
useful  experiment  to  not  only  validate  the  design,  but  also  to  overcome  the  difficulties  of 
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building  a  structure  that  can  withstand  an  atmosphere  of  external  pressure  without 
leaking. 

A  study  of  the  diffusivity  of  various  membranes  would  be  necessary  to  detennine 
a  suitable  skin  for  the  geodesic  structure.  A  vacuum  LTA  vehicle  could  potentially 
require  less  servicing  than  a  helium  LTA  vehicle  because  the  membrane  preventing  the 
equalization  of  pressure  would  require  pores  smaller  than  the  critical  diameter  of  an 
Oxygen  molecule  of  2.8  A  as  opposed  to  the  critical  diameter  of  a  Helium  molecule  of  2 
A.  [1].  The  pressure  gradient  across  the  membrane  is  a  factor  in  diffusivity  rates  as  well, 
so  that  the  relative  diffusivity  rates  between  a  traditional  helium  filled  membrane  which 
has  very  little  pressure  gradient,  and  the  theoretical  air  displacing  vacuum  membrane 
which  has  an  atmosphere  pressure  gradient  are  not  clear  without  further  study. 

An  optimization  should  be  performed  on  the  skin  where  W/B  is  minimized  by 
varying  skin  thickness  where  W/B  is  affected  by  volume  loss  due  to  membrane 
displacement  and  weight  due  to  thickness.  The  skin  will  most  likely  need  to  be  composed 
of  a  low  diffusivity  membrane  such  as  Mylar  and  a  high  strength  reinforcement  cloth 
such  as  Zylon. 

The  rotating  cylinder  concept  could  be  further  studied  by  researching  the  torque 
and  acrit  for  gap  sizes  between  1  and  0.  Depending  on  the  results  of  this  study  rotating 
cylinder  LTA  vehicle  filled  with  helium  could  be  evaluated  for  feasibility  using  von 
Kannan  efficiency. 
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Appendix  A:  Matlab  Code 


A.l  Optimization  Routine  (LTA  opt.m) 

%LTA  fmincon  and  Genetic  Algorithm 

%Written  by  1st  Lt  Justin  Mason,  Capt  Trent  Metlen,  and  Brian  Cranston 
clear;  close  all;  clc; 

~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k 

%  This  program  is  for  MECH  620  optimization  of  a  lighter  than  air  (LTA) 
%  icosahedron  structure. 

O, 

o 

■ k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

tic 

%  Scaling  factors  for  design  variables  for  optimization  routine: [f  r  t] 
scaling  =  [1  1  1 ] ; 

%%  Optimization  This  section  picks  which  optimization  routine  to  run 

%  Set  upper  and  lower  bounds 
ub=[l.l  le-1  le-2 ]. /scaling; 
lb=[0.9  le-5  le-6] . /scaling; 

%  Set  initial  guess 

%  Note:  this  is  only  used  in  fmincon 
xo  =  [1  4e-3  2e-4] . /scaling; 

optimization  =  'fmincon'; 

%optimization  =  'genetics- 

switch  optimization 
case  'genetic' 

%%  Genetic  Algorithm  Optimization 

%  Run  an  optimization  routine  with  a  genetic  algorithm 

%  A  =  [f  r  t] 

%  g3  =  inequality  const 

A=  [ 0  -1  .5];  b=0 ; 
nvars  =  3; 

IntCon  =  1; 

options=gaoptimset (' display ' , ' iter ' , ' Populations ize ; , 10 ) ; 

[x, fval, exitflag, output, population, scores]  =  ga(@(x) 

Cost (x,  scaling) , nvars , A, b,  [ ] ,  [ ] , lb, ub, @ (x) 

Constraints (x, scaling, krieging) , options) ; 
x  =  x . *scaling; 

case  ' fmincon ' 

%%  fmincon  Optimization 

%  This  section  runs  a  fmincon  optimization  routine 
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%  A  =  [f  r  t] 

%  g3  =  inequality  const 

A=  [ 0  -1  .5];  b= [ 0 ] ; 

options=optimset ( ' algorithm '  , 'interior-point', 'display', 'iter'); 

[x,fval]  =  fmincon (@ (x)  cost (x, scaling) , xo, A, b,  [],[],  lb,  ub,  @ (x) 

Constraints (x, scaling, krieging) , options) ; 

x  =  x . *scaling; 

end 

toe 


A.2  Optimization  Sub-Routine:  Cost  Function  (Cost.m) 

function  [  cost  ]  =  Cost (  x, scaling) 

load  geometryl .mat 
x=x . *scaling; 
f  =  round (x ( 1 ) )  ; 
rcb  =  x  ( 2 )  ; 
teb  =  x ( 3 ) ; 

%  Get  sphere  parameters 
[par]  =  Sphere_parameters (x) ; 

%  Constants 
rhocb  =  par  ( 1 ) ; 
rhosk  =  par (2 ) ; 
rs  =  par  ( 3 ) ; 
tsk  =  par ( 4 ) ; 

%  leb  =  par ( 5 ) ; 

%  Volume=par ( 6 ) ; 

rhoair=l . 225;  %kg/mA3,  density  of  air 
buoyancy=rhoai r* Volume; 

cost  =  (rhocb* ( 60*pi*tcb*rcb* (leb) *fA2) +Area*tsk*rhosk) /buoyancy; 
end 


A.3  Optimization  Sub-Routine:  Sphere  Parameters  (Sphere  parameters.m) 

function  [  par  ]  =  Sphere_parameters (  x  ) 

%  Capt  Trent  Metlen,  1st  Lt  Justin  Mason,  and  Brian  Cranston 


f  =  round (x ( 1 ) ) ; 
r  =  x  ( 2  )  ; 
t  =  x  ( 3 )  ; 

load  geometryl .mat 

%average  lengths  of  geodesic  spheres  f=l:10  for  rs=0.33 

store_lcb  =  [0.3470  0.1926  0.1313  0.0993  0.0797  0.0666 

0.0572  0.0501  0.0445  0.0401] ; 

Icb  =  pi* ( (r+ . 5*t) A4- (r- . 5*t) A4 ) /4 ;  %area  moment  of  inertia 

rhosk  =  1400;  %kg/mA3  Ref  1 
rs=0 . 33 ; 
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tsk=l . 6e-5; 

lcb  =  store_lcb (round (f) ) ; 
sigcbcy  =  le9;  %UHM  carbon  epoxy 

Ecb  =  450e9;  %modulus  of  elaticity  UHM  carbon  epoxy 
rhocb  =  1522;  %density  UHM  carbon  epoxy 

par  =  [rhocb; rhosk; rs ; tsk; lcb; Volume; sigcbcy; Ecb; lcb] ; 
end 

A.4  Optimization  Sub-Routine:  Constraints  (Constraints.m) 

function  [  G, ceq  ]  =  Constraints (  x, scaling) 

x  =  x.*scaling; 
f  =  round (x ( 1 ) ) ; 
rcb  =  x  ( 2 ) ; 
tcb  =  x  ( 3 )  ; 

[par]  =  Sphere_parameters (x) ; 
sigcbcy  =  par (7); 
rs  =  par ( 3 ) ; 

gl=  abs (sigmax) -sigcbcy; 

[sigmax  g2]  =  Analyze  beam ( f, rcb, tcb) ; 
g4=  rcb+ . 5*tcb-sqrt (3) *rs/4/f ; 

%  nonlinear  constraints  gl,  g2,  g4 
G  =  [gl;  g2  ;  g4]  ; 

ceq  =  [  ] ; 
end 


A.5  Optimization  Sub-Routine:  Finite  Element  Program 

function  [  sigmax  g2  ]  =  Analyze  beam(f,r,t) 

%  IstLt  Justin  Mason  and  Dr  Black's  Finite  Element  Class  2012 
%  Create  and  analyze  FE  models  from  *.dat  files 

%%  Create  the  *.dat  file 
Datcreate ( [round (f)  r  t] ) 

%%  Read  *.dat  file 
model=Datread_beam ( ' icos . dat ' ) ; 

%  disp('Read  in  dat  file'); 

%%  Initialize 
do  f  =  6  ; 

N  =  model . ngrid*dof; 

K  =  zeros (N, N) ; 

%%  Generate  global  K  matrix 

for  element  index  =  1  :  model. ncbar  ; 

EID  =  model . cbar . EID (element  index); 
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PID  =  model . cbar . PID (EID) ; 

MID  =  model .pbar .MID (PID) ; 

G1  =  model . cbar . Gl (EID)  ; 

G2  =  model . cbar . G2 (EID) ; 

A  =  model . pbar .A ( PID) ; 

E  =  model .matl . E (MID) ; 

L  =  sqrt ( (model . grid . XI (G2 ) -model . grid . XI (Gl ))  A2+ ..  . 
(model . grid . X2 (G2 ) -model . grid .X2 (Gl) ) A2+. . . 
(model . grid . X3 (G2 ) -model . grid .X3 (Gl) ) A2)  ; 


I^zz=model .pbar. II (PID)  ; 
I  yy=model . pbar . 12 ( PID)  ; 
J=model . pbar . J (PID)  ; 
nu=model .matl . NU (MID)  ; 
G=E/ (2* (1+nu) ) ; 


X  = 
Y_1 
Y_2 
Y_3 
Y_4 
Z_1 
Z_2 
Z_3 
Z_4 
S  = 
elemK 


A*E/L; 

=  12*E*I_zz/LA3; 
=  6*E*I_zz/LA2; 

=  4*E*I_zz/L; 

=  2*E*I_zz/L; 

=  12*E*I_yy/LA3; 
=  6*E* I_yy/LA2 ; 

=  4*E*I_yy/L; 

=  2*E*I_yy/L; 

G* J/L; 


X 

0 

0 

0 

0 

0 

-X 

0 

0 

0 

0 

0 

0 

Y  1 

0 

0 

0 

Y  2 

0 

-Y  1 

0 

0 

0 

Y  2 

0 

0 

Z  1 

0 

-Z  2 

0 

0 

0 

-Z  1 

0 

-Z  2 

0 

0 

0 

0 

S 

0 

0 

0 

0 

0 

-S 

0 

0 

0 

0 

-Z  2 

0 

Z  3 

0 

0 

0 

Z  2 

0 

Z  4 

0 

0 

Y  2 

0 

0 

0 

Y  3 

0 

-Y  2 

0 

0 

0 

Y  4 

X 

0 

0 

0 

0 

0 

X 

0 

0 

0 

0 

0 

0 

-Y  1 

0 

0 

0 

-Y  2 

0 

Y  1 

0 

0 

0 

-Y  2 

0 

0 

-Z  1 

0 

Z  2 

0 

0 

0 

Z  1 

0 

Z  2 

0 

0 

0 

0 

-s 

0 

0 

0 

0 

0 

S 

0 

0 

0 

0 

-Z  2 

0 

Z  4 

0 

0 

0 

Z  2 

0 

Z  3 

0 

0 

Y  2 

0 

0 

0 

Y  4 

0 

-Y  2 

0 

0 

0 

Y  3 

O, 

o 


o, 

o 


%  Create  and  apply  transformation  matrix  T 


T=zeros (12,12)  ; 

V— [model . cbar .XI (EID) ; model . cbar . X2 (EID) ; model . cbar . X3 (EID) ] ; 
Vx= [model . grid .XI (G2 ) -model . grid . XI (Gl )  ,  .  .  . 
model . grid . X2 (G2 ) -model . grid . X2 (Gl )  ,  .  .  . 
model . grid . X3 (G2 ) -model . grid . X3 (Gl )  ]  ; 

Vz=cross (Vx, V) ; 

Vy=cross (Vz , Vx)  ; 

Vx  e=Vx/norm (Vx)  ; 

Vy_e=Vy/norm (Vy)  ; 

Vz_e=Vz/norm(Vz)  ; 

delta  =  [Vx_e; Vy_e; Vz_e]  ; 

z=zeros (3,3) ; 
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T= [delta  z  z  z;  z  delta  z  z;z  z  delta  z;z  z  z  delta]; 

TelemK=T ' *elemK*T; 

O,  _ 

o 

%  Place  elemK  into  global  K 

o,  _  _  _ _ _  _ 

o 

Gl_index=f ind (model . grid . ID==G1 )  ; 

G2  index=f ind (model . grid . ID==G2 ) ; 
dof sl= (Gl_index-1 ) *dof+l : G1  index* dof; 
dof s2= (G2_index-1 ) *dof +1 : G2  index* dof; 

K (dof si , dof si ) =K (dof si , dof si ) +TelemK ( 1 : dof , 1 : dof )  ; 
K(dofs2,dofsl) =K (dof s2 , dof si ) +TelemK (dof +1 : 2* dof , 1 : dof ) ; 

K (dof si , dof s2 ) =K (dof si , dof s2 ) +TelemK ( 1 : dof , dof +1 : 2  *dof ) ; 
K(dofs2,dofs2) =K (dof s2 , dof s2 ) +TelemK (dof +1 : 2*dof , dof +1 : 2*dof ) 


end 

%  disp (' Generated  global  K  matrix'); 

%%  Create  force  vector 
F  =  zeros (N,l);  %Initialize 

for  ctr  =  1  :  model. nforce  %Cycle  through  all  F0RCE1  cards 
G_index=f ind (model . grid . ID==double (model . force . G (ctr ) ) ) ; 
fdof s= (G_index-1 ) *dof +1 : (G_index-1 ) *dof +3 ; 

Vxf = [model . grid . XI (model . force . G2 (ctr) ) - 
model . grid . XI (model . force . G1 (ctr)  )  ;  .  .  . 

model . grid . X2 (model . force . G2 (ctr) ) - 
model . grid . X2 (model . force . G1 (ctr)  )  ;  .  .  . 

model . grid . X3 (model . force . G2 (ctr) ) - 
model . grid . X3 (model . force . G1 (ctr)  )  ]  ; 

Vx  e  force=Vxf /norm (Vxf )  ; 

%Force  vec=model . force . F 

F (fdofs) =F (fdofs) +Vx  e_f orce*model . force . F (ctr) ; 
end 

%  disp (' Created  Force  vector'); 

%%  GRDSET  Cards 
%  cdof= [ ] ; 

cdof=false (N, 1) ; 

for  ctr  =  1  :  model . ngrdset  %Cycle  through  all  grdset  cards  if 
model . ngrdset 

for  ctr2=l: length (model . grdset . PS ) 

dim=str2double (model . grdset . PS (ctr2 ) )  ; 
cdof ( (0 : (model .ngrid-1) ) *dof+dim) =true; 

%  cdof= [cdof, ( (0 : (model .ngrid-1) ) *dof+dim) ] ; 

end 

end 

%%  SPC  Cards 

for  ctr  =  1  :  model. nspc  %Cycle  through  all  SPC1  cards 
for  dim_ctr=l : length (model . spc . C { ctr } ) 

dim=str2double (model . spc . C { ctr } (dim_ctr ) ) ; 

G_idx=f ind (model . grid . ID==double (model . spc . G1 (ctr) ) )  ; 
cdof ( (G_idx-1 ) *dof+dim) =true; 

end 
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end 

udof=~cdof  ; 


%%  Compute  Displacements 

D  =  zeros (N,l);  %create  displacement  vector 

D(udof)  =  K (udof ,  udof ) \F (udof )  ;  %  solve  for  displacements 

F(cdof)  =  K(cdof, :)*D;  %  compute  reaction  forces 

%  disp  (  '  '  ) 

%  disp('  Displacements') 

%  disp ('Node  T1  T2  T3  R1  R2 

R3  '  )  ; 

%  disp  (  ' -  -  -  -  -  - 


for  ctr  =  1  :  model. ngrid 

row  =  (ctr-1) *dof+l  :  ctr*dof; 

%  disp (sprintf ( ' %4d  %10.2e  %10.2e  %10.2e  %10.2e  %10.2e 

%10 . 2e ' , . . . 

%  full (model .grid. ID (ctr) ) , D (row, 1) ) ) ; 

end 


%  %%  Reaction  Forces 
%  disp  (  '  '  ) 

%  disp('  Reaction  Forces') 

%  disp ('Node  FI  F2  F3  Ml  M3 

M3  '  )  ; 

%  disp  (  ' -  -  -  -  -  - 


o, 

o 

%  for  ctr  =  1  :  model. ngrid 
%  row  =  (ctr-1) *dof+l  :  ctr*dof; 

%  disp ( sprintf (' %4d  %10.2e  %10.2e  %10.2e  %10.2e  %10.2e 

%10 . 2e ' ,  . . . 

%  full (model .grid. ID (ctr) ) , F (row, 1) ) ) ; 

%  end 

%%  Stress  and  strains 
%  disp  (  '  '  ) 

%  disp  (  '  Stress  '  ) 

%  disp (' Element  SigmaX  Taut  T  Tau  Y  Tau 

z' )  ; 

%  disp  ( ' -  -  -  -  - 

' )  ; 

sig  all  =  zeros ( 1 , length (model . ncbar) ); g2  =  sig  all; 
for  ctr=l :model . ncbar 

EID  =  model . cbar . EID (ctr) ; 

PID  =  model . cbar . PID (EID)  ; 

MID  =  model .pbar. MID (PID)  ; 

G1  =  model . cbar . G1 (EID) ; 

G2  =  model . cbar . G2 (EID)  ; 

A  =  model . pbar .A ( PID) ; 

E  =  model .matl . E (MID)  ; 

L  =  sqrt ( (model . grid . XI (G2 ) -model . grid . XI (G1 )) A2+ .. . 
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(model . grid .  X2 (G2 ) -model . grid .X2 (Gl) ) A2+. . . 

(model . grid . X3 (G2 ) -model . grid .X3 (Gl) ) A2) ; 

I  zz=model .pbar . II (PID)  ; 

I  yy=model . pbar . 12 ( PID)  ; 

J=model . pbar . J (PID)  ; 
nu=model .matl . NU (MID) ; 

G=E/ (2* ( 1+nu) ) ; 

T=zeros  ( 12 , 12 )  ; 

V= [model . cbar .XI (EID) ; model . cbar . X2 (EID) ; model . cbar . X3 (EID) ]  ; 

Vx= [model . grid .XI (G2 ) -model . grid . XI (Gl )  ,  .  .  . 
model . grid . X2 (G2 ) -model . grid . X2 (Gl )  ,  . . . 
model . grid . X3 (G2 ) -model . grid . X3 (Gl ) ]  ; 

Vz=cross (Vx, V) ; 

Vy=cross (Vz , Vx) ; 

Vx  e=Vx/norm (Vx)  ; 

Vy  e=Vy/norm (Vy) ; 

Vz  e=Vz/norm (Vz ) ; 

delta  =  [Vxe; Vy_e; Vz_e] ; 

z=zeros (3,3) ; 

T= [delta  z  z  z;  z  delta  z  z;z  z  delta  z;z  z  z  delta]; 

TelemK=T ' *elemK*T; 

Gl  index=f ind (model . grid . ID==G1 )  ; 

G2  index=f ind (model . grid . ID==G2 )  ; 
dofsl= (Gl_index-1) *dof+l : Gl_index*dof ; 
dof s2= (G2_index- 1 ) *dof +1 : G2_index*dof ; 

D_locall= [D (dof si ) ; D (dof s2 ) ]  ; 

D  local=T*D  locall; 

ul=D_local ( 1 , 1 ) ;  vl=D_local (2 , 1 ) ;  wl=D_local(3,l); 
u2=D_local ( 7 , 1 ) ;  v2=D_local ( 8 , 1 ) ;  w2=D_local ( 9 , 1 ) ; 

theta_xl=D_local (4,1) ;  theta_yl=D_local (5, 1) ;  theta_zl=D_local (6,1) 
theta_x2=D_local (10, 1) ;  theta_y2=D_local (11, 1) ; 
theta  z2=D  local (11,1); 
x=0  ; 
y=r  ; 
z=0  ; 

cy=2 ; cz=2 ; 

M_z=E*I_zz* ( (-6/LA2+12*x/LA3) *vl+ (-4/L+6*x/LA2) *theta_zl+ ( 6/LA2- 
12*x/LA3) *v2+ (-2/L+6*x/LA2) *theta_z2) ; 

V_z=E*I_zz* ( 12/LA3*vl+6/LA2*theta_zl-12/LA3*v2+6/LA2*theta_z2 ) ; 
M_y=E*I_yy* ( (-6/LA2+12*x/LA3) *wl- (-4/L+6*x/LA2) *theta_yl+ ( 6/LA2- 
12*x/LA3) *w2- (-2/L+6*x/LA2) *theta_y2) ; 

V_y=E*I_yy* (12/LA3*wl-6/LA2*theta_yl-12/LA3*w2-6/LA2*theta_y2) ; 
Normal=A*E/L* (u2-ul+v2-vl+w2-wl ) ; 

Torsion=G*J* (theta  x2-theta  xl+theta_y2-theta_yl+theta  z2- 
theta_zl ) /L; 

Sigma  x=Normal/A-M  z*y/I  zz-M  y*z/I_yy; 

Tau_t=Torsion* (r+t/2) /J; 

Tau_y=cy*V_y/A; 

Tau  z=cz*V  z/A; 
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%  disp ( sprintf ( ' %4d  %10.2e  %10.2e  %10.2e 

%10 . 2e '  ,  .  .  . 

%  full (EID) , full (Sigma_x) , full (Tau_t) , full (Tau_y) , full (Tau_z) ) ) 

sig  all(2*ctr-l)  =  Sigma  x/2-sqrt ( (Sigma  x/2 ) A2+Tau_t A2 ) ; 

%  disp ( sprintf (' sigl  %10 . 2e ' , sigall (2*ctr-l) ) ) 

g(2*ctr-l,l)  =  Sigma_x*A-piA2*E*I_zz/L/ . 5; 
x=L; 

M_z=E*I_zz* ( (-6/LA2+12*x/LA3) *vl+ (-4/L+6*x/LA2) *theta_zl+ ( 6/LA2- 
12*x/LA3) *v2  + (-2/L+6*x/LA2) *theta_z2) ; 

V_z=E*I_zz* (12/LA3*vl  +  6/LA2*theta_zl-12/LA3*v2  +  6/LA2*theta_z2)  ; 
M_y=E* I_yy* ( (-6/LA2+12*x/LA3) *wl- (-4/L+6*x/LA2) *theta_yl+ ( 6/LA2- 
12*x/LA3) *w2- ( -2 /L+6*x/LA2 ) *theta_y2) ; 

V_y=E*I_yy* (12/LA3*wl-6/LA2*theta_yl-12/LA3*w2-6/LA2*theta_y2) ; 
Normal=A*E/L* (u2-ul+v2-vl+w2-wl ) ; 

Torsion=G*J* (theta  x2-theta  xl+theta_y2-theta_yl+theta  z2- 
theta_zl )  /In¬ 
signia  x=Normal/A-M  z*y/I  zz-M  y*z/I_yy; 

Tau_t=Torsion* (r+t/2) /J; 

Tau_y=cy*V_y/A; 

Tau  z=cz*V_z/A; 

%  disp ( sprintf ( '  %10.2e  %10.2e  %10.2e 

%10 . 2e ' ,  .  .  . 

%  full (Sigma  x) , full (Tau  t) , full (Tau  y) , full  (Tau  z) ) ) ; 

g(2*ctr,l)  =  -Sigma_x*A-4*piA2*E*I  zz/LA2; 

sig  all (2*ctr)  =  Sigma  x/2-sqrt ( (Sigma  x/2)A2+Tau  tA2); 

%  disp ( sprintf (' sigl  %10 . 2e ' , sigall (2*ctr) ) ) 

end 

g2  =  max ( full (g) ) ; 
sigmax  =  min (sig  all); 
sigmean=mean ( sig  all) ; 

%  disp ('Done  analyze_beam ' ) 
end 

A.6  Optimization  Sub-Routine:  Datcreate  (Datcreate.m) 

function  [  ]  =  Datcreate (  x  ) 

%  IstLt  Justin  Mason 

%  This  program  writes  a  .dat  file  for  optimization  routine  of  an 
%  icosahedron 

2-9-  ★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★ 

O  O 

%  THIS  FUNCTION  MUST  BE  RUN  WITHIN  LTAgeod.mat 
%  inputs:  rcb,  tcb,  XYZ,  beam,  &  force 
%%  Define  Constants 

[XYZ,  beam,  force,  nV,  centroid]  =  LTAgeod2 (x) ; 

f  =  x  ( 1 )  ; 
rcb  =  x  ( 2 ) ; 
tcb  =  x  ( 3 )  ; 
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%area  moment  of  inertia 


A  =  2*pi*rcb*tcb;  %area 

11  =  pi* ( (rcbt . 5*tcb) A4- (rcb- . 5*tcb) A4 ) /4 ; 

12  =  II;  %area  moment  of  inertia 

J  =  pi* ( (rcbt . 5*tcb) A4- (rcb- . 5*tcb) A4 ) /2 ;  %torsional  constant 
E  =  450e9;  %modulus  of  elasticity  of  carbon  epoxy 

nu  =  0.33;  %poisson's  ratio  of  carbon  epoxy 

rho  =  1522;  %density  of  carbon  epoxy 

%%  Open  file  and  create  header 
fileName  =  'icos.dat'; 
fid  =  fopen (fileName, ' w' ) ; 

header  =  { ' ID  ICOS  OPT  '  ;  ' SOL  1 0 1 CEND TITLE  =  ICOS  OPT'; 
'SUBTIT  =  MECH  62 O' ;' LABEL  =  ICOSAHEDRON SPC  =  1'; 

'LOAD  =  1 ' ;  ' D I S P  =  ALL' ;' STRESS  =  ALL' ;' BEGIN  BULK';}; 
for  ctr  =  1 : length (header ) 

fprintf (fid, s treat (header { ctr } ,  \n ' ) ) ; 

end 


%%  Write  XYZ  Locations 

fprintf ( fid,  ' $  Geometry\n ' ) ; 

for  ctr  =  1 :  length (XYZ ) 

fprintf (fid,  ' GRID, %d, ,%.4f,%.4f,%.4f\n',XYZ(ctr,  :) )  ; 

end 


%%  Write  beam  Conectivity 

fprintf ( fid,  ' $  Beam  Connectivity\n ' ) ; 

for  ctr  =  1 : length (beam) 

idx  =  [ f ind (XYZ ( : , 1 ) ==beam (ctr , 2 ) )  f ind (XYZ ( : , 1 ) ==beam (ctr, 3) ) ] ; 
v  =  XYZ (idx, 2:4); 
orient  =  null (v) ; 

fprintf (fid,  ' CBAR, %d, 1 , %d, %d, % . 4f , % . 4f , % . 4f \n ' , beam (ctr ,  :)  ,  orient) 

end 

%%  Write  beam  and  Material  Properties 

fprintf ( fid, ' $  Beam  and  Material  Properties\n ' ) ; 

fprintf (fid, ' PBAR, l,l,%.10e,%.10e,%.10e,%.10e,0.0\n',A,Il,I2,J); 

fprintf (fid,  MAT1 , 1 , % . 4f , ,%.4f, % . 4f \n '  ,  E,  nu,  rho)  ; 

%%  Write  Boundary  Conditions 

fprintf ( fid, ' $  Boundary  Conditions\n ' ) ; 

fprintf (fid, 'GRDSET, ,,,,,, \n' ) ; 
fprintf (fid, ' SPC 1,1, 123456, %d\n ' , nV) ; 
for  ctr  =  nV+1 : nV+length (centroid) +1 

fprintf (fid,  'SPC 1,  1,  123456,  %d\n '  ,  XYZ (ctr ,  1 )  )  ; 

end 

for  ctr  =  1 : length (force) 

fprintf (fid,  ' FORCE 1 ,  l,%d,%.4f,%d,%d\n'  ,  force (ctr ,  : ) ) ; 
end 
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%  fprintf (fid, ' FORCE 1, 1, 1,1000, 1, 2\n ' ) ; 
fprintf (fid, ' ENDDATA ' ) ; 

fclose (fid) ; 

end 


A.7  Optimization  Sub-Routine:  Geodesic  Sphere  Geometry  (LTAgeod2.m) 

function  [  XYZ,  beam,  force,  nV,  centroid]  =  LTAgeod2 (  x  ) 

%Capt  Trent  Metlen  10/10/2012 

%  Frequency  of  triangles  per  major  triangle  (1  for  icosahedron) 

f =x ( 1 ) ; 
rcb  =  x  ( 2 ) ; 
tcb  =  x  ( 3 )  ; 
r=0 . 33 ; 

%  Angles  in  spherical  coordinates  used  to  determine  icosahedron 
vertices 

phil=0. 4636476090008; 
phi2=pi/ 2-phil ; 
theta l=2*pi/5; 

%  Atmospheric  Conditions 

Pa=l . 5*101325;  %pascals,  sea  level  pressure  (safety  factor  1.5) 

%  Icosahedron  verticies  in  spherical  coordinates 

Vp=[0  pi/2;0  phil;thetal  phil ;  2*thetal  phil ; 3*thetal  phil ; 4*thetal 
phil;  0.5*thetal  -phil ; 1 . 5*thetal  -phil ; 2 . 5*thetal  -phil ; 3 . 5*thetal  - 
phil;  4.5*thetal  -phil;0  -pi/2]; 

%  Basic  information  on  geodesic  shape 
nV=10*fA2+2;  %number  of  vertices 
nt=20*fA2;  %number  of  triangles 
ne=30*fA2;  %number  of  edges 

%  Find  all  vertices  of  first  5  major  triangles  for  geoesic  sphere 
V=  [r*ones (nV, 1 )  ones(nV, 2)]; 

V  ( 1 , 2 : 3 ) =Vp ( 1 , :); 

V ( nV ,2:3) =Vp (12,  :)  ; 

K=1  ; 

for  1=1 :f 

for  J=1:I*5 

V( (K+J) , 2 : 3) = [ (J-l) *2*pi/ (1*5)  pi/2-I*phi2/f ] ; 

end 

J; 

K=K+ J; 

end 
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%  Find  all  vertices  of  next  10  major  triangles 
for  1=1 :f 

for  J=l:5*f 

V(  (K+J)  , 2 : 3)  =  [I* (pi/ (5*f) )  +  (2*pi/ (5*f )  )  * (J-l)  phil-2*phil*I/f ] 

end 

J; 

K=K+ J; 

end 

%  Find  all  vertices  of  last  5  major  triangles 
for  1=1 : f 
if  I<f 

for  J=1 : 5* (f-I) 

V(  (K+J) , 2 :3)  =  [pi/5+ (2*pi/ (5* (f-I) ))  * (J-l)  -pi/2+phi2-phi2*I/f  ] ; 
end 
K=K+ J; 

end 

end 

%  Convert  from  spherical  coordinates  to  cartesian  coordinates 

[x  y  z] =sph2cart (V ( : , 2) ' , V ( : , 3) ' , V ( : , 1) ' ) ; 

leg=sqrt ( (x (1) -x (2) )A2+(y(l)-y(2) )A2+(z (1) — z (2) )A2) ; 

%  Create  xyz  matrix  in  cartesian  coordinates 
XYZ= [x '  y'  z']; 

%  Find  beam  matrix  for  icosahedron  (f=l) 
if  f==l 

beam=[l  2;1  3;1  4;1  5;1  6; 

2,11;2,7;3,7;3,8;4,8;4,9;5,9;5,10;6,10;6,11; 

7, 12; 8, 12; 9, 12; 10, 12; 11, 12; 

2,3;3,4;4,5;5,6;6,2; 

7, 8;  8, 9; 9, 10; 10, 11; 11, 7]  ; 
n=length (beam) ; 
beam=[(l:n)'  beam] ; 

else  %for  more  complex  geodesic  spheres  ( f > 1 ) 

%  Split  up  variables  into  matrices  of  planes  in  sphere 
p=2*f+f+l ; 

plane  =  cell (1, 2*f+l) ; 
plane] 1}=[1  XYZ(1,:)]; 

if  f>l 

cntr=l ; 
for  I  =  2 : f 

cntr2= (1-1) *5; 

plane { I } = [ (cntr+1 : cntr2+cntr )  '  (XYZ (cntr+1 : cntr2+cntr ,  : ) ) ] 
cntr=cntr2+cntr ; 

end 

end 

for  I=f+1 : 2*f+l 
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plane { I } = [ (cntr+1 : cntr+f *5 )  '  XYZ (cntr+1 : cntr+f *5 ,  : ) ]  ; 
cntr=cntr+f*5; 

end 
J=0  ; 

for  I=p-  (f-1)  :p-l 
J=J+1 ; 

plane { I } = [ (cntr+1 : cntr+ ( f- J) *5 ) '  XYZ (cntr+1 : cntr+ ( f- J) *5 , : ) ] ; 
cntr=cntr+ (f-J) *5; 

end 

plane{p}= [cntr+1  XYZ (end, :)]; 

%  Determine  number  of  planes 
[spl, sp2] =size (plane) ; 

%  Create  connections  between  top  point  and  plane  2 
cntr=l : cntr ; 

beam=[l  2;1  3;1  4;1  5;1  6];  %beam= [counter (n)  nodel  node2] 
n=6;  %n  is  the  beam  counter 

%  Create  connections  between  coplaner  nodes  (creates  horizontal 
rings ) 

for  1=2 : sp2-l 

[si, s2]=size (plane {1,1}) ; 
for  J=l:sl-1 

beam(n, :)=[plane{l,I} (J,l)  plane{l,I} (J+1,1) ] ; 
n=n+l ; 

end 

beam(n,  :)  =  [plane{l,I}  ( si , 1 )  plane{l,I}  (1,1)]; 
n=n+l ; 

end 

%  Create  connections  between  planes  in  triangular  cap  starting  plane 
2  to  3 

e=[l  1  1];  ledge  point  connections  (points  lying  on  edge  of  major 
triangles ) 

c=[l  1];  %center  point  connections  (points  not  lying  on  edges  of 
major  triangles) 

for  1=2 :f 

%  B=matrix  of  e's  and  c's  where  e ' s=l  and  c's=0 

B=[l  zeros (1,1-2)  1  zeros(l,I-2)  1  zeros(l,I-2)  1  zeros(l,I-2) 

1  zeros (1,1-2)  ]  ; 

[si,  s2] =size (B) ; 

%  A=matrix  of  points  where  the  first  row  are  start  points,  second  row 
%  end  points,  will  be  added  to  beam  matrix  when  loop  is  complete 

%  first  point  connections 

A=  [plane { 1 , 1 }  ( 1 , 1 ) *e;  plane { 1 , 1  +  1 }  (end, 1 )  plane { 1 , 1  +  1 }( 1 , 1 ) 
plane { 1 , 1  +  1 }  (2, 1)  ]  ; 
cntr=2 ; 
for  J=2:s2 
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if  B ( J) >0 
m=e; 

else 

m=c; 

end 

%  All  other  point  connections  on  plane 
A=[A(1,:)  plane { 1 , I } ( J, 1 ) *m;  A(2,:) 

(plane{l,I+l} (cntr: cntr+length (m) -1 , 1 ) ) ' ] ; 
cntr=cntr+length (m) -1; 

end 

[s4,  s3] =size (A) ; 

%  Beam  matrix  is  updated  by  A  matrix 

beam (n : n+s3-l ,  : )  =  [A ( 1 ,  : )  '  A (2 ,  : )  ' ]  ; 
n=n+s3 ; 

end 

%  Create  connections  from  bottom  of  triangular  cap  to  bottom  of 
middle  triangles 

m= [ 1  1]; 

for  1=1 : f 
s2=f*5; 

A= [plane{l,f+I} ( 1 , 1 ) *m;  plane{l,f+I+l} (end, 1 ) 
plane {l,f+I+l}(l,l)]; 
cntr=l ; 
for  J=2:s2 

A=[A(1,:)  plane { 1 ,  f+I }  (J, 1) *m;  A(2,:) 

(plane{l,f+I  +  l}  (cntr: cntr+1 , 1 ) )  ' ]  ; 
cntr=cntr+l ; 

end 

[s4,  s3] =size (A) ; 

beam(n:n+s3-l,  : )  =  [A (1 ,  :)  '  A(2,  :)  '  ]  ; 
n=n+s3 ; 

end 

%  Create  connections  from  bottom  of  middle  triangles  to  last  plane 
%  before  endpoint 

e= [ 1 ] ; 
c= [1  1]; 

for  I=l:f-1 

B=[l  zeros(l,f-I)  1  zeros(l,f-I)  1  zeros(l,f-I)  1  zeros(l,f-I) 
1  zeros ( 1 ,  f-I )  ]  ; 

[si, s2] =size (B) ; 

A= [plane{l,2*f+I} (end, 1 ) *c;  plane {l,2*f+I+l} (end, 1 ) 
plane {l,2*f+I+l}  (1,1) ]  ; 
cntr=l ; 
for  J=l:s2-1 
if  B ( J) >0 
m=e; 

else 

m=c; 

end 
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A= [A (1 , : )  plane {l,2*f+I} (J,l)*m;  A(2,:) 

(plane {l,2*f+I+l} (cntr: cntr+length (m) -1 , 1 ) ) ' ] ; 
cntr=cntr+length (m) -1 ; 

end 

[s4,  s3] =size (A) ; 

beam (n : n+s3-l ,  : )  =  [A ( 1 ,  :)  '  A  (2 ,  :)  ']; 
n=n+s3 ; 

end 

%  Create  endpoint  connections 

A=[  (planed,  3*f}  (1:5,1)  )  (plane  { 1 , 3*f+l }( 1 , 1  ))*[  1  111  1]]; 

beam(n:n+4, : ) = [A (1 , :) '  A(2, :) ']; 
n=n+4 ; 

beam=[(l:n)'  beam]; 

end  %ends  if  loop  which  was  necessary  in  the  case  of  f=l 

%  show  3D  plot  of  geodesic  sphere 
%  Plot  shape 
%  tic 

%  figure ( ) ; 

%  for  1=1 : n 

%  plot3 (XYZ (beam ( 1 , 2 : 3 ) ,1) , XYZ (beam ( 1 , 2 : 3 ) ,2) , XYZ (beam ( 1 , 2 : 3 )  ,3) ) 

o, 

o 

%  hold  on 

%  end 
%  toe 

%  Find  Triangles  Making  up  Geodesic  Sphere 
triangle=[0  0  0]; 
for  1=1 : nV 

ind  =  f ind (beam ( : , 2 ) ==I )  ; 

ind  =  [ind; f ind (beam ( : , 3) ==I ) ] ; 

A=beam ( ind, : ) ; 
points=unique (A ( : , 2 : 3 ) ) ; 

points ( find (points==I ))=[] ;  %remove  point  of  interest  from  list  of 
points 

for  J=1 : length (points ) 

ind  =  f ind (beam ( : , 2 ) ==points ( J) ) ; 

ind  =  [ind; find (beam ( : , 3) ==points ( J) ) ] ; 

A=beam ( ind, : ) ; 

points 2 =unique (A ( : , 2 : 3 ) ) ; 

points2 ( find (point s2==I ) )  =  []  ; 

points2 (find (point s2==points ( J) ) ) = [ ] ; 

B=zeros (size (points) )  ; 
for  K=1 : length (points2  ) 

Bl=points==points2 (K) ; 

B=B+B1 ; 

end 

point  vert ices=non zeros (B . *points ) ; 
triangle= [triangle; I* (ones (size (point_vertices ) ) ) 
points ( J) * (ones (size (point_vertices ) ) )  point_vertices ] ; 
end 

end 
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triangle  ( 1 ,  :)  =  []; 

for  1=1 : length (triangle) 

triangle (I, : ) =sort (triangle (I, : ) ) ; 

end 

triangle=unique (triangle, ' rows ' ) ; 
for  1=1 : length (triangle) 

centroid ( I ,  : ) =mean ( [ (XYZ (triangle (1,1),:));  (XYZ (triangle (1,2),  : ) )  ; 
(XYZ (triangle (1,3),:))]); 

%calc  height  of  each  tetrahedron 
height ( I ) =norm (centroid ( I ,  : ) )  ; 

%calc  area  of  each  triangular  face 
areat (I) =0 . 5*norm(cross ( (XYZ (triangle (1,3),:)- 
XYZ (triangle (1,1),:)),  (XYZ (triangle (1,3),:) -XYZ (triangle (1,2),  :)))); 
end 

%calculate  volume  of  each  irregular  tetrahedron 
volumet= (1/3) *areat . *height; 

%find  volume  and  surface  are  of  total  structure 
Volume=sum (volumet)  ; 

Area=sum (areat) ; 

%  center  point  index  in  xyz  matrix 
cp=nt+nV+l ; 

%calc  force  vector  force= [point  mag  ptl  pt2] 

force=  [0  0  0  0 ] ; 

for  1=1 : length (triangle) 

%  find  vectors  for  each  triangle  to  determine  area 

vl2=0 . 5* (XYZ (triangle (1,2), : ) -XYZ (triangle (1,1),:)); 

vlc=centroid (I,  : ) -XYZ (triangle (1,1),  :  )  ; 

vl3=0 . 5* (XYZ (triangle (I,  3)  ,  : ) -XYZ (triangle (1,1),:)); 

v21=0 . 5* (XYZ (triangle (1,1), : ) -XYZ (triangle (1,2),:)); 

v2c=centroid ( I ,  : ) -XYZ (triangle (1,2),  :  )  ; 

v23=0 . 5* (XYZ (triangle (1,3), : ) -XYZ (triangle (1,2),:)); 

v31=0 . 5* (XYZ (triangle (1,1), : ) -XYZ (triangle (1,3),:)); 

v3c=centroid (I, : ) -XYZ (triangle (1,3), :); 

v32=0 . 5* (XYZ (triangle (1,2), : ) -XYZ (triangle (1,3),:)); 

%  Calc  force  magnitude=  Pressure*area  where  area  is  1/2  cross  product 
%  of  vectors  (triangle) .  Two  triangles  per  irregular  quadrilateral. 

%  Three  quadrilatersals  make  up  the  total  area  of  a  triangular  face 
%  with  the  pressure  from  each  quadrilateral  acting  as  an  equivalent 
%  force  at  the  closest  vertex  of  the  triangular  face.  Each 
%  quadrilateral  is  defined  by  four  points,  one  at  the  vertex  of  the 
%  triangular  face,  one  at  the  centroid,  and  the  remaining  two  halfway 
%  along  each  edge  of  the  triangular  face  that  meets  the  vertex  at  which 
%  the  force  will  act. 

magn=Pa* . 5* [norm (cross (vl2 , vie) ) +norm (cross (vie, vl3 ) ) ; 
norm (cross (v21 , v2c) ) +norm (cross (v2c,  v2  3 ) )  ; 
norm (cross (v31 , v3c) ) +norm (cross (v3c, v32 ) ) ] ; 

%  magnitude  = [point  mag  ptl  pt2] 

%  [point  of  force  application,  magnitude  of  force, 

%  starting  point  of  force  direction,  ending  point  of  force 
direction] 
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%  start  point  is  centroid  of  triangular  face,  end  point  is  origin 
force= [force;  triangle ( I , 1 )  magn(l)  nV+I  cp; 
triangle ( I , 2 )  magn(2)  nV+I  cp; 
triangle ( I , 3 )  magn(3)  nV+I  cp]  ; 

end 

force  (1,  :)  =  []; 

%  Find  average  length  of  beams 
lcbi=zeros (1, length (beam) ) ; 
for  1=1 : length (beam) 

lcbi ( I ) =norm (XYZ (beam (1,3),  : ) -XYZ (beam (1,2),  : ) )  ; 

end 

lcb=mean (lcbi) ; 

save  geometryl .mat  Volume  Area  lcb; 

%  Add  centroids  and  origin  to  xyz  matrix 
XYZ= [XYZ; centroid; 0  0  0]; 

%  Add  index  to  xyz  matrix 
XYZ= [ ( 1 : cp) '  XYZ] ; 

End 


A.8  Optimization  Sub-Routine:  .dat  File  Reader  (Datread  beam.m) 


function  [model]  =  Datread  beam (filename) 

%  IstLt  Justin  Mason  and  Dr  Black's  2012  Finite  Element  Class 
%  Simple  *.dat  file  reader  with  very  limited  capabilities 
%  Warning  -  use  at  your  own  risk! 

%%  Open  file 

fid  =  fopen (filename) ;  %Open  dat  file 


%%  Initialize  Variables 


model . ncrod=0 ; 

O, 

o 

Number 

of 

rod  element 

cards 

model . ncbar=0  ; 

o, 

o 

Number 

of 

bar  element 

cards 

model . ngrid=0  ; 

Q. 

O 

Number 

of 

grid  point  cards 

model . nprod=0  ; 

O, 

O 

Number 

of 

rod  property 

cards 

model . npbar=0 ; 

O, 

O 

Number 

of 

bar  property 

cards 

model . nmat=0 ; 

O, 

o 

Number 

of 

material  cards 

model . nspc=0 ; 

o, 

o 

Number 

of 

single  point 

constraint 

model . nf orce=0 ; 

o, 

o 

Number 

of 

force  cards 

model . ngrdset=0 ; 

o, 

o 

Number 

of 

gridset  cards 

model . grid . ID  =  sparse (zeros)  ; 
model . grid . CP  =  sparse (zeros) ; 
model . grid . XI  =  sparse (zeros) ; 
model . grid . X2  =  sparse (zeros) ; 
model . grid . X3  =  sparse (zeros) ; 
model . grid . dofl  =  sparse (zeros) ; 
model . grid . dofn  =  sparse (zeros) ; 
model . crod . EID  =  sparse (zeros)  ; 
model . crod . PID  =  sparse (zeros) ; 
model . crod . G1  =  sparse (zeros) ; 
model . crod . G2  =  sparse (zeros) ; 

model . cbar . EID  =  sparse ( zeros ); %Beam  Elements 


cards 
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model . cbar . PID  =  sparse (zeros) ; 
model . cbar . G1  =  sparse (zeros) ; 
model . cbar . G2  =  sparse (zeros) ; 
model . cbar . XI  =  sparse (zeros) ; 
model . cbar . X2  =  sparse (zeros) ; 
model . cbar . X3  =  sparse (zeros) ; 
model . prod . PID  =  sparse (zeros) ; 
model . prod . MID  =  sparse (zeros)  ; 
model. prod. A  =  sparse (zeros) ; 

model . pbar . PID  =  sparse ( zeros ); %Beam  Mat ' 1  Properties 

model . pbar . MID  =  sparse (zeros) ; 

model. pbar. A  =  sparse (zeros) ; 

model.pbar.il  =  sparse (zeros) ; 

model . pbar . 12  =  sparse (zeros) ; 

model. pbar. J  =  sparse (zeros) ; 

model . pbar . NSM  =  sparse (zeros) ; 

model .matl . MID  =  sparse (zeros) ; 

model. matl.E  =  sparse (zeros) ; 

model. matl. G  =  sparse (zeros) ; 

model .matl . NU  =  sparse (zeros) ; 

dof s=6 ; 

%%  Read  Executive  Control 

read_more  =  true;  %  intialize  logical  var  to  stop  reading  dat  file 
first  =  textscan ( f id,  ' %s 1 delimiter '  ');  %read  line  up  to  first 

space 

first  string  =  char ( first { 1 , 1 }) ;  %  Transform  cell  array  to  character 
string 

while  (read  more)  %  While  logical  read  more  is  true  keep  reading 
switch  first  string  %  Switch  on  the  first  string  that  is  read  in 

case  {'ID'}  %Read  in  the  ID  for  the  job 

cell  =  textscan ( fid,  '%[A\n]',l);  %read  the  rest  of  the  line 
model. ID  =  char (cell { 1 ,  1 })  ;  %convert  to  a  character  string 

case  {'SOL'}  %Read  in  the  solution  sequence  to  be  executed 
cell  =  textscan ( fid,  ’%u',l);  %read  the  rest  of  the  line 
model. sol  =  uint32 (cell { 1 , 1 } )  ;  %convert  to  an  integer 

case  { ' CEND ' }  %Exit  when  done  with  executive  control 
read  more  =  false;  %A11  done  reading 
break 


end 

first  =  textscan ( fid,  ' %s ', 1 ,' delimiter ', '  ');  %read  line  up  to  the 

first  space 

first_string  =  char ( first { 1 ,  1 })  ;  %  Transform  cell  array  to  character 

string 
end 

%  disp (' Completed  reading  Executive  Control');  %Dump  to  the  screen  an 
update 

%%  Read  Case  Control 


148 


read  more  =  true;  %  intialize  logical  var  to  stop  reading  dat  file 
first  =  textscan ( f id,  ' %s 1 , 1 delimiter '  ');  %read  line  up  to  the 

first  space 

first  string  =  char (first { 1 , 1 }) ;  %  Transform  cell  array  to  character 

string 

while  (read  more)  %  While  logical  read  more  is  true  keep  reading 
switch  first  string  %  Switch  on  the  first  string  that  is  read  in 

case  {'TITLE'}  %Read  in  the  title 

cell  =  textscan ( fid,  '%[A\n]',l);  %read  the  rest  of  the  line 
temp  =  char (cell { 1 , 1 }) ;  %convert  to  a  character  string 
model. title  =  temp (2 : length (temp) ) ;  %trim  off  the  '=' 

case  {'SUBTIT'}  %Read  in  subtitle 

cell  =  textscan ( fid,  '%[A\n]',l);  %read  the  rest  of  the  line 
temp  =  char (cell { 1 , 1 }) ;  %convert  to  a  character  string 
model . subtitle  =  temp (2 : length (temp) ) ;  %trim  off  the  '=' 

case  {'LABEL'}  %Read  in  the  label 

cell  =  textscan ( fid,  '%[A\n]',l);  %read  the  rest  of  the  line 
temp  =  char (cell { 1 , 1 }) ;  %convert  to  a  character  string 
model. label  =  temp (2 : length (temp) ) ;  %trim  off  the  '=' 

case  {'SPC'}  %Read  in  the  SPC  identifier 

cell  =  textscan ( fid,  ’%s%u',l);  %read  the  rest  of  the  line 
model. spc_ID  =  uint32 (cell { 1 , 2 } )  ;  %convert  to  an  integer 

case  {'LOAD'}  %Read  in  the  Load  identifier 

cell  =  textscan ( fid,  *%s%u',l);  %read  the  rest  of  the  line 
model . load_ID  =  uint32 (cell { 1 , 2 } )  ;  %convert  to  an  integer 

case  {'DISP'}  %Read  in  the  Disp 

cell  =  textscan ( fid,  '%[A\n]',l);  %read  the  rest  of  the  line 
temp  =  char (cell { 1 , 1 }) ;  %convert  to  a  character  string 
model. disp  =  temp (2 : length (temp) ) ;  %trim  off  the  '=' 

case  {'STRESS'}  %Read  in  the  Stress 

cell  =  textscan ( fid,  '%[A\n] ',1);  %read  the  rest  of  the  line 
temp  =  char (cell { 1 , 1 }) ;  %convert  to  a  character  string 
model. stress  =  temp (2 : length (temp) ) ;  %trim  off  the  '=' 

case  {'BEGIN'}  %Exit  when  done  with  case  control 
read  more  =  false;  %A11  done  reading 
break 


end 

first  =  textscan ( fid,  ' %s ', 1 ,' delimiter ', '  ');  %  Read  line  up  to  the 

first  space 

first_string  =  char ( first { 1 , 1 }) ;  %  Transform  cell  array  to  character 

string 
end 

%  disp (' Completed  reading  Case  Control');  %Dump  to  the  screen  an 
update 


149 


%%  Read  Bulk  Data 

read  more  =  true;  %  intialize  logical  var  to  stop  reading  dat  file 
first  =  textscan ( f id,  ' %s 1 , 1 ,' delimiter ',',') ;  %read  line  up  to  the 
first  comma 

first  string  =  char ( first { 1 , 1 }) ;  %  Transform  cell  array  to  character 

string 

while  (read  more)  %  While  logical  read  more  is  true  keep  reading 
switch  first  string  %  Switch  on  the  first  string  that  is  read  in 


case  {'GRID1, 'EGRID'}  %Read  in  a  grid  card 

model. ngrid  =  model. ngrid  +  1;  %Increment  the  number  of  grid 


points 

cell  =  textscan  ( fid,  '  %• 
' emptyValue '  ,  0 )  ; 
GN  =  uint32 (cell { 1 , 1 } ) ; 
model . grid . ID (model . ngr 
model . grid . CP (GN) 
coordinate  system 

model . grid . XI (GN) 
coordinate 

model . grid . X2 (GN) 
coordinate 

model . grid . X3 (GN) 
coordinate 

model . grid . dofl (GN) 
model . grid . dof n (GN) 


%u%f%f%f ' , 1, 'delimiter' , 

%read  the  rest  of  the  line 
%Store  grid  ID 
d)  =  GN; 

=  uint32 (cell { 1 , 2  } )  ;  %Store 

=  double (cell { 1 , 3  })  ;  %Store 

=  double (cell { 1 , 4  })  ;  %Store 

=  double (cell { 1 , 5  }) ;  %Store 

=  (model . ngrid-1 ) *dofs+l ; 

=  (model . ngrid) *dofs; 


1 

2 

3 


case  {'CROD'} 

model. ncrod  =  model. ncrod  +  1; 

cell  =  textscan (fid,  '%u%u%u%u’,l,' delimiter 

'emptyValue',  0);  %read  the  rest  of  the  line 
EID  =  uint32 (cell { 1 , 1 } ) ; 

model . crod . EID (model . ncrod)  =  EID;  %Store  element  ID 
model . crod . PID (EID)  =  uint32 (cell { 1 , 2 } ) ;  %Store  the 

property  ID 

model . crod . G1 (EID)  =  uint32  (cell { 1 , 3  } ) ;  %Store  grid  ID 

1 

model . crod . G2 (EID)  =  uint32 (cell { 1 , 4 } ) ;  %Store  grid  ID 

2 


case  {'PROD'} 

model. nprod  =  model. nprod  +  1; 

cell  =  textscan (fid,  ' %u%u%f ' , 1 , ' delimiter 

'emptyValue',  0);  %read  the  rest  of  the  line 
PID  =  uint32  (cell { 1, 1 } )  ; 

model . prod . PID (model . nprod)  =  PID;  %Store  property  ID 
model . prod . MID ( PID)  =  uint32 (cell { 1 , 2 } ) ;  %Store  material  ID 
model .prod. A (PID)  =  double (cell { 1 , 3 }) ;  %Store  cross  sectional 

area 

case  { ' CBAR ' } 

model. ncbar  =  model. ncbar  +  1; 

cell  =  textscan ( fid,  ' %u%u%u%u%f %f %f ', 1 ,' delimiter ' 
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r  •  • 


' emptyValue ' ,  0);  %read  the  rest  of  the  line 
EID  =  uint32 (cell { 1, 1 } ) ; 

model . cbar . EID (model . ncbar)  =  EID;  %Store  element  ID 
model . cbar . PID (EID) 
property  ID 

model . cbar . G1 (EID) 

1 


uint32  (cell { 1 , 2 } ] 
uint32  (cell { 1 , 3  }  ] 


%Store  the 
%Store  grid  ID 


2 


model . cbar . G2 (EID) 


model . cbar .XI (EID) 
component  xl  of  orientation  vector 
model . cbar . X2 (EID) 
component  x2  of  orientation  vector 
model . cbar . X3 (EID) 
component  x3  of  orientation  vector 


uint32 (cell{l,4} 
double (cell { 1 , 5 } 
double (cell { 1 , 6 } 
double (cell { 1 , 7 } 


;  %Store  grid  ID 
;  %Store 
;  %Store 
;  %Store 


case  { ' PBAR ' } 

model. npbar  =  model. npbar  +  1; 

cell  =  text scan (fid,  ' %u%u%f%f%f%f%f',l, ; delimiter 
'emptyValue',  0);  %read  the  rest  of  the  line 
PID  =  uint32 (cell { 1, 1 } )  ; 

model . pbar . PID (model . npbar )  =  PID;  %Store  property  ID 
model . pbar . MID ( PID)  =  uint32 (cell { 1 , 2 } )  ;  %Store  material  ID 
model . pbar .A ( PID)  =  double (cell { 1 , 3 }) ;  %Store  cross  sectional 


area 

model .pbar. II (PID) 

inertia 

model . pbar . 12 (PID) 

inertia 

model . pbar . J (PID) 

inertia 


=  double (cell { 1 , 4  }) ; 
=  double (cell { 1 , 5  }) ; 
=  double (cell {1,6}); 


%Store  1st  moment  of 
%Store  2st  moment  of 
Store  polar  moment  of 


model . pbar . NSM ( PID)  =  double (cell { 1 , 7 }) ;  %Store  non- 

structural  mass 


case  { ' MAT1 ' } 


model. nmat  =  model. nmat  +  1; 

cell  =  text scan (fid,  '%u%f%f%f',l, ' delimiter 

'emptyValue',  0);  %read  the  rest  of  the  line 
MID  =  uint32 (cell { 1 , 1 } ) ; 

model .matl .MID (model .nmat)  =  MID;  %Store  material  ID 
model .matl . E (MID)  =  double (cell { 1 , 2 }) ;  %Store 

modulus 

model .matl . G (MID)  =  double (cell { 1 , 3 }) ;  %Store 

modulus 


model .matl . NU (MID) 

ratio 


=  double (cell { 1 , 4 }) ;  %Store 


Youngs 

Shear 

Poissons 


case  {'SPC1'} 

model. nspc  =  model. nspc  +  1;  %Increment  the  number  of  SPC  cards 
cell  =  text scan (fid,  ' %u%s%u ' , 1 , ’ delimiter 

'emptyValue',  0);  %read  the  rest  of  the  line 
model . spc . SID (model . nspc)  =  uint32 (cell { 1 , 1 } ) ;  %Store  SPC  ID 
model . spc . G1 (model . nspc)  =  uint32 (cell { 1 , 3 } ) ;  %Store  node  to  be 
constrained 
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Store 


model . spc . C {model . nspc }  =  char (cell { 1 , 2 } ) ; 

constrained  components 


case  { ' F0RCE1 ' } 

model. nforce  =  model. nforce  +  1;  %Increment  the  number  of 


F0RCE1  cards 

cell  =  text scan (fid,  ' %u%u%f %u%u ' , 1, ' delimiter 


'emptyValue',  0);  %read 
model . force . SID (model . nforce) 

the  rest  of  the  line 
=  uint32 (cell { 1 , 1 } ) ; 

%Store 

forcel  ID 

model . force . G (model . nforce) 

=  uint32 (cell { 1 , 2 } ) ; 

%Store 

ID 

model . force . F (model . nforce) 

=  double (cell { 1 , 3  })  ; 

%Store 

magnitude  of  force 

model . force . G1 (model . nforce) 

=  uint32 (cell { 1 , 4  } )  ; 

%Store 

begining  node  ID  for  vector 

model . force . G2 (model . nforce) 

=  uint32 (cell { 1 , 5  } ) ; 

%Store 

ending  node  for  ID 


node 


case  { ' GRDSET ' } 

model . ngrdset  =  model . ngrdset  +  1;  %Increment  the  number  of 
GRDSET  cards 

cell  =  text scan (fid,  '%u%u%u%u%u%u%s',l, ' delimiter 
' emptyValue ' ,  0);  %read  the  rest  of  the  line 
model . grdset . CP  =  uint32 (cell { 1 , 2 } ) ;  %Coordinate  system  for 
grid  points 

model . grdset . CD  =  uint32 (cell { 1 , 6 } ) ;  %Coordinate  system  for 
displacements 

model.grdset.PS  =  char (cell {1,7}); 


case  'ENDDATA'  %Exit  when  done  with  data 
read  more  =  false;  %A11  done  reading 
break 


end 

first  =  textscan ( f id,  ' %s ', 1 ,' delimiter ;  %read  line  up  to  the 
first  comma 

first  string  =  char ( first { 1 , 1 }) ;  %  Transform  cell  array  to  character 

string 
end 

%  disp (' Completed  reading  Bulk  Data');  %Dump  to  the  screen  an  update 

fclose(fid);  %Open  dat  file 

end 

A.9  Von  Karman  Efficiency  for  Rotating  Cylinder  (LTAheliumcylinder.m) 

%  LTAheliumcylinder 

%  Capt  Trent  Metlen 

clc;  close  all;  clear  all; 

global  P  rho  ro  mfact  LIE  Vs  Ms  sm  d  J 

%*************************atm0sphere*********************************** 
rho=0.889;  %  kg/mA3,  density  of  air  at  10k  feet 
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P=1500*70;  %Pa,  Pressure  and  1.5  safety  factor 
mu=l . 68*10 A-5;  %dynamic  viscosity  of  air 

^-k  k  ~k  k  k  k  ~k  k  k  k  ~k  k  -k  k  k  k  k  k  k  k  -k  k  k  k  k  k  ~k  k  k  k  ~k  k  k  k  ~k  k  k  k  ~k  k  k  k  ~k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  ~k  k  k  k  ~k  k  k 

v=.3;  %  Poisson's  ratio 

sm=5 . 67*10A9;  %Pa  tensile  strength  of  Hextow  lm7 

E=2.76*10All;%2.76*10All;%modulus  of  Hextow  lm7  *10A11 

Ec=l . 2 *  1 0 A 1 0 ;  %transverse  modulus  elasticity  of  UHM  carbon  fiber 

ss=128*10A6;  %shear  strength  of  AS4  hextow  fiber 

dhe=0.056;  %kg/mA3  density  of  helium  at  10k  ft 

da=rho-dhe; %density  of  air  -  density  of  helium 

d=1522 . 3945; %density  of  AS4/3501-6 

Vc=0  ; 

mc=0  ; 

K=l; 

syms  cdi 

options  =  optimset ( 1  Display of f 1 ) ;  %  Turn  off  display 

options2=optimset ( 1  Display' ,  ' off  '  ,  'Algorithm' ,  ' trust-region- 

reflective ')  ;  %  Turn  off  display 

mub=.0011;  %bearing  friction  coefficient  Ref  5 

db= . 15; %bearing  inner  diameter  (m)  Ref  6 

to= .0045; 

for  J=1 : 20 

ro=8;%m,  radius 

u=10*J;%  m/s, tip  speed  of  cylinder 

Re=rho*u . *ro/mu;  %reynolds  number  for  cylinder  in  free  flow 
for  counter=l : length (Re) 

Rel=Re (counter) ; 

F=solve ( 10 A ( ( 1/ (cdi A . 5 ) + . 6 ) /4 . 07 ) -Rel*cdi A  .  5 ) ;  %find  cd 
cd (counter) =double (F) ;  %convert  symbolic  to  numeric 

end 

L=ro* [1 : 1 : 10] ; %m,  length  of  cylinder 
As=4*pi*roA2 ; %mA2 ,  area  of  spherical  endcaps 
for  1=1 : length (L)  %find  required  thickness  and  mass  for  each  L  and 
R 

eta=ro/to; 

Kc=0 .4233+79. 9779*etaA-l-12759 . 6621*etaA-2+755633 . 25*etaA-3; 
z=L (I) A2*sqrt (l-vA2) / (ro*to) ; 

Ac ( I ) =L ( I ) *2*pi* (ro) ; %mA2 ,  surface  area  of  cylinder 
M=0 . 5*rho*ro*Ac ( I ) *cd* (u . A2 ) ;  %  Nm,  moment 

t (I) = (M*zA (1/4) / (2*pi*ro*Kc*E) ) A . 5; %m,  required  thickness  of 
%cylinder  shell  due  to  torque 
to=t (I )  ; 

tl=M/ (2*pi*roA2*ss) ;  %thickness  required  due  to  shear  stress 
if  tl>to  %check  local  buckling  condition,  is  ultimate  strength 
%greater  than  local  stress? 

t  (I) =tl; 

end 

ts (I) =M. / (2*pi* ( ( . 15+ro) / 2)  A2*ss) ; %m,  thickness  of 
spherical  endcaps 

Ms ( I ) =As*ts ( I ) *d; %kg,  mass  of  sphere 
Vs (I) =4/3*pi*roA3; %mA3,  volume  of  sphere 
me (I) =d*Ac (I) *t (I) +Ms (I) ; %kg,  mass  of  cylinder  with  spherical 

%endcaps 

Vc (I) =L (I) *pi* (ro) A2+Vs (I) ; %mA3,  volume  of  cylinder  with 
%spherical  endcaps 
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Lr ( I ) =L ( I ) / (ro) ; %Length/radius  of  cylinder 

end 

wc=mc*9.8;  %N,  weight  of  spherical  endcaps 

wac=Vc*da* 9 . 8 ;  %N,  buoyant  force 

Lift=wac-wc;  %N,  lift  force  of  cylinder 

fb= . 5 . * (Lift) ;  %N,  bearing  load  weight  of  1/2  payload 

Cm=  .  14  6  .  /  (Re  . A .  2  )  ;  %Moment  Coefficient 

Pb=0 . 5*mub*db* (u . /ro) . *fb;  %watts,  power  required  due  to 
bearing  friction 
w=u/ro; 

Ps=l /2*Cm. *rho . *u . A3 . *ro . A2 ;  %watts,  power  required  to  rotate 
disks  (ends  of  cylinders) 

Pr=w.*M  %watts,  power  required  to  rotate  cylinder 
Px=0 . 75*Pr+Ps+Pb;  Ref  1  75%  power  requirement 
U=u/2.5;  %m/s,  forward  velocity  omegacrit=2 . 5  gap  size=2*r 
evk=Px . / (U*Lift) ;  %von  Karman  efficiency 

[Evk ( J) , ind] =min (evk (find (evk>0) ) ) ;  %find  minimum  evk  assuming 
no  engine  weight  (so  that  minimization  happens  even  if  vehicle  is 
infeasible) 

[~, ind] =find (evk==Evk ( J) ) 

Acm ( J) =Ac (ind) ; 

rm ( J) =Lr (ind) ;  %length/radius  at  optimal  evk 
tm ( J) =t (ind) ;  %thickness  at  optimal  evk 
rom(J) =ro;  %radius 

wci ( J) =wc (ind)  ;  %weight  of  cylinder  at  optimal  evk 
waci ( J) =wac (ind) ;  %buoyancy  of  cylinder  at  optimal  evk 
omega=w;  %rotational  velocity  required  for  equilibrium 
Vtip(J)=u;  %tip  speed  of  cylinder  in  m/s 
power (J) =Px (ind) ;  %power  required  to  spin  cylinder  NM/ s 
Ux ( J) =U; 

end 

volc=2*pi* (roA3) .  *rm;  %(mA3)  volume  of  two  cylinders  used  to  compare 
to  conventional  airship 

Lif t=waci-wci ;  %lbs  buoyancy-weight 
Lm=wci . /waci; 

zeta=43 . 15*10A6;  %(Nm/kg)  energy  density  of  fuel 
%calculate  C27  evk 

fuel=4000;  %kg  of  fuel  used  by  C-27J  to  go  4074  km 
dist=4459616;  %distance  traveled  (m) 
velocity=150 ;  %m/s  cruise  speed  of  C-27J 
c2 7=fuel*zeta* velocity/ dist; 
ec27=c27/ (10000*9. 8* velocity)  ; 

x=0 :max (rom) ; 

rc=6;  %ft  radius  of  cylindrical  gondola  cabin 
lc=60;  %ft  length  of  gondola  cabin 

vc=pi*rcA2*lc*0 . 0283168466  ;  %volume  of  gondola  mA3 
mc=ll*vc;  %mass  of  gondola  (kg)  based  on  11  kg/mA3  Ref  11 
mb=26*4;  %mass  of  bearings  (kg)  Ref  A  6 

figure  ( ) 
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lambda=2.5;  %fineness  ratio 
Vol=volc (5 ) ; 

r= ( (3/4) *Vol/ (pi* lambda) ) A (1/3) ; 
rhol0=.406;  %(kg/mA3)  Ref  2 

Dc= . 5*rho*Ux . A2 . *pi*rcA2*0 . 012 ;  %Drag  of  cabin  Cd=.012  Ref  12 
b=Vol* . 85*da*9 . 8;  % (N)  buoyancy  volume  of  lifting  gas  is  .85  of  total 
airship  volume 

CD=.0169  %CD  of  BLC  airship  Ref  3  .0162  +.0007  for  tail  surfaces 
U= (1/1 . 9438) * [10 : 10 : 350] ;  %m/s,  velocity 

k=l/.855;  %.59  is  propeller  efficiency  with  45%  improvement 
Ap=pi* (r/4 . 64) A2;  %  mA2,  propeller  area  (propeller  radius=  max 
radius/4 . 64 ) 
range=1852000;  %  m 

zeta=43 . 15*10A6;  %Nm/kg,  energy  density  of  Jet  A 

Px= (k*Dc . A2 . / (2*rho*Ap*Ux) ) ;  %Power  required  to  overcome  cabin  drag 
Lx=2*Lift- (2*1. 5*9. 8* (power+Px) ) .*0.00134- 
2* 9 . 8*range . /Ux . * (power+Px) . *3/ zeta  ;  %N,  Useful  Load  of  rotating 
cylinder  LTA  vehicle  1.5  kg/hp  Ref  6 
ex=3*power./ (Ux.*Lx);  %Von  Karman  Efficiency  of  rotating  cylinder  LTA 
%vehicle  assuming  33  percent  efficient  engines 
for  1=1 : length (U) 

t (I) — range/U (I) ;  %time  enroute 
D (I) =. 5*rholO*U (I) A2*VolA (2/3) *CD;  % (N)  Drag 
P ( I ) =k*sqrt (D ( I ) A3/ (2*rho*Ap) ) ;  %Power  required 
L (I) =b* .7- (1 .7*9. 8*P (I) ) *0 .00134-9. 8*t (I) *P (I) *3/zeta  ;  % (N) 
useful  load  assuming  30%  structure  to  useful  load  ratio  and  1.7  kg/hp 
engine  and  ducted  fan  weight  as  well  as  33%  engine  efficiency  and 
energy  density  of  Jet  A 

e  (I) =3*P (I) / (U (I) *L (I) ) ;  %Von  Karman  efficiency  at  each  velocity 
assuming  33  percent  efficient  engines 
end 

%von  Karman  efficiency  of  zeppelin  NT 
eznt=147000/ (1900*9 . 8*34)  ; 

%Gabrielli  von  Karman  limit 
A=  .000175; 

V=1 : 10 : 1010;  %mph 
ef f k=A*V; 

%  Plot  von  Karman  efficiency  graph 

loglog (V*0 . 869, effk, U* 1 . 944, e, Ux*l . 94  4 , ex,  ' k ' , velocity* 1 .944,ec27,  ' b  * ' , 
34*1 . 944, eznt, ' r .' )  %mph* . 8 6 9=knots  m/s*l . 944=knots 
text (velocity*1.6,ec27*1.5, ' C-27 ' ) 
hold  on 

data3= [  1/3.9  22.4;  1/.42  75.8;  1/.21  8.9;  1/.4  15.6;  1/.62  7.5; 
1/36.68  26.8;  1/1.86  29.5;  1/2.54  290.5;  1/320  8.9;  1/21  24.6];%Ref  1 
ef f=data3 (:,1) ./(data3(:,2)*2.237)  ; 

scatter(data3(:,2)*1.944,data3(:,l), ' . ')  %m/s*l ,944=knots 
text ( 1 . l*Ux ( 10 ) *1 . 944 , ex ( 10 ) , 'Rotating  Cyl ' ) 

text ( 1 . 1*0 . 8*U ( length (U))*1.944,0.8*e( length (U) ) , ' BLC  Airship ' ) 

text (34*1.6,eznt*1.5, ' Zeppelin  NT ' ) 

text (data3 (1,2) /1.5,data3 (1,1) , 'Motorcycle' ) 

text (data3 (2,2) /1.2,data3 (2,1)  ,  ' N99VE ' ) 

text (data3 (3, 2) /3, data3 (3, 1) , 'Automobile  in  City') 

155 


text (data3(4,2)/3,data3(4,l),  ' Commuter  Train ' ) 
text (data3 (5,2) ,data3 (5,1) , 'Bus' ) 
text (data3 ( 6, 2 ) /2 , data3 ( 6, 1 )  , 'Freight  Train') 

text (data3 (7,2) /2,data3 (7,1) , 'Airship  1936' ) 
text (data3 (8,2) *1 . 6, data3 (8, 1) *1.5,  '747  '  ) 
text (data3(9,2)/1.5  , data3 ( 9 , 1 ) , 'Tanker') 
text (data3 ( 10 , 2 ) , data3 ( 10 , 1 )  ,  ' Truck ' ) 

xlabel ({ 'Velocity  (knots)  ' },  ' font size ' ,14) ; 
ylabel({'Von  Karman  Ef f iciency ' } , ' fontsize ' , 14 ) ; 
legend('Von  Karman  Limit', ' BLC  Airship', 'Rotating  Cylinder  LTA 
Vehicle ' , ' Location ' , ' Best ' ) 

axis ( [1  1000  10 A-4  10]) 

A.10  Optimization  of  Isogrid  Stiffened  Sphere  (LTAsphere.m) 

%LTAsphere 
%Capt  Metlen 

clc;  clear  all ;  close  all; 
global  t  r  v  sc  E  P  rs 
%  atmospheric  properties 

da=1.225;  %kg/mA3,  density  of  air  at  sea  level 
P=1500*101 . 325  %Pa,  air  pressure  sea  level  and 

%  properties  of  skin 
%  Dupont  E-130-x 

E=.5*895e9;  %Pa,  Young's  modulus  of  elasticity 
v=0.3;  %  Poisson's  ratio 
sc=1786e6;  %Pa,  compressive  strength 
s=3445e6;  %Pa,  tensile  strength 
dens=1522;  %kg/mA3,  density 

r=0.33;  %m,  radius  of  sphere 

ti=linspace (. 000001 , 0 . 0001 , 100 ) ;  %m,  thickness  of  skin 
xo=[0.013  0.0001  0.0004];  %initial  guess  xo=[d  b  h] 
for  1=1 : length  (ti) 
t=ti  (I)  ; 

x(I,:)  =  fsolve (@LTAspheref ,  xo,  optimset (' Display ',  ' of f ')) ;  %solve  for 
[d  b  h]  at  each  thickness 
xo— x ( I , : ) ; 
end 

d=abs (x  ( : , 1) )  '  ; 
b=abs (x  ( : , 2) )  '  ; 
h=abs (x ( : , 3) )  '  ; 

a=2*h/sqrt (3) ;  %m,  triangle  leg  length 
te=ti+3*d . *b . /h;  %m,  smeared  thickness 
tes=ti+d . *b . /h; 

scrit=P*r./ (2*te);  %stress  in  cross  section  due  to  uniform  external 
pressure 
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te=te-te . * ( scrit>sc) + ( scrit>sc) . *P*r . / (2*sc) ;  %make  sure  material 
compressive  strength  not  exceeded.  If  stress  is  greater  than  material 
compressive  strength  then  thickness  is  defined  by  material  compressive 
strength,  otherwise  it  is  defined  by  buckling  of  sphere 

alpha=b. *d. / (ti . *h) ; 
del=d .  /ti; 

beta= (3*alpha .  * (ones (size (del) ) +del) . A2+ (ones (size (alpha) ) +alpha) . * (one 
s (size (alpha) ) +alpha . *del . A2 ) ) .  A  .  5 ; 

ts=ti . *beta ./ (ones (size (alpha) ) +alpha) ;  %  m,  equivalent  thickness 
Es=E* ( (ones (size (alpha) ) +alpha) . A2) . /beta;  %  Pa,  equivalent  stiffness 

tsmooth=r*sqrt (P*sqrt (l-vA2) / ( . 8*E) )  ;  %  m,  thickness  required  of 
unstiffened  sphere 

scrits=P*r . / (2*tsmooth)  ;  %stress  in  unstiffened  sphere 

tsmooth=tsmooth- tsmooth* (scrits>sc) + (scrits>sc) . *P*r./ (2*sc)  %make  sure 
material  compressive  strength  not  exceeded 
A=4 *pi . *r A2 ; %  (mA2)  area  of  sphere 

V=4 /3*pi* ( r+ . 5*te)  . A3 ;  %  mA3,  volume  of  air  displaced  by  sphere 

W=A*te*dens* 9 . 8 ;  %  N,  Weight  of  skin 

bo=9.8*V*da;  %  N,  buoyant  force 

Mb=W./bo;  %weight  to  buoyancy  ratio 

L= (bo-W) ;  % (N)  Lift 

Wsmooth=A*tsmooth*dens* 9 . 8 ;  %  N,  Weight  of  skin 

Mbsmooth= (Wsmooth . /bo) ;  %weight  to  buoyancy  ratio  of  unstiffened  sphere 
beam=min (find (a . /d>5) )  ; 

figure ( ) ; 

plot (ti, Mb, ti (beam: end) , Mb (beam: end) , ' r ' , ti, Mbsmooth) ; 
xlabel('Skin  Thickness,  meters');  ylabel ( ' W/B ' ) 

%title (' Weight  to  Buoyancy  Ratio  of  Be  Blade  Stiffened  Sphere') 
legend('Grid  Stiffened  Sphere',  'Unstiffened  Sphere',  'Location',  'Best') 
%axis([ti(l)  ti (length (ti) )  min (Mb) *.8  max (Mb)]) 

function  F  =  LTAspheref (x) 
global  t  r  v  sc  E  P 
d=abs(x(l));  %blade  depth 
b=abs(x(2));  %blade  width 
h=abs(x(3));  %triangle  height 
a=2*h/sqrt (3) ;  %triangle  leg  length 
te=t+3*d*b/h;  %smeared  thickness 
tes=t+d*b/h; 

scrit=P*r/ (2*te);  ^stress  in  cross  section  due  to  uniform  external 
pressure 

alpha=b*d/ (t*h) ; 
del=d/t; 

beta= (3*alpha* (1+del) A2+ (1+alpha) * ( l+alpha*del A2 ) ) A . 5; 
ts=t*beta/ (1+alpha) ; 

Es=E* ( (1+alpha) A2) /beta; 

if  (a/d)>5  %  solve  buckling  of  beam  for  rib  buckling 
F= [scrit- (piA2) *E* (bA2) / (12*aA2) ; P-0 . 8*Es*tsA2/ ( (rA2) *sqrt (1- 
vA2 ) ) ; scrit-5*  (piA2) *E* ( (t/a) A2) / (12* (l-vA2) ) ] ; 
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Else  %  solve  buckling  of  plate  for  rib  buckling 

F=[scrit-0.456* (pi A2 ) *E* (bA2 ) / ( 12* (dA2 ) * (l-vA2) ) ; P- 
0 . 8*Es* (tsA2) / ( (rA2) *sqrt (l-vA2) ) ;  scrit-5* (piA2)*E*((t/a)A2)/(12*(l- 
vA2 ) ) ] ; 
end 

A.ll  Partial  Vacuum  Study  (LTAspherepv.m) 

%Capt  Metlen 
%LTAspherepv 
clc;  clear;  close  all; 

ni=logspace (0, 3, 1000)  ; 

E=0 . 0795*450*10A9;  %0.0795  from  LTAsphere  gives  E*  at  optimal  stiffener 
setting 

patm=l . 5*1 . 01325*10 A5;  %1  bar  =  10A5  Pa  (N/m2) 
mu=.3;  %Poisson's  ratio 
r=l;  %m  radius  of  sphere 

rhos=1522;  %kg/mA3  density  of  surface  of  sphere 
rhoair=l . 225;  %kg/mA3  density  of  air  at  sea  level 
for  1=1 : length (ni ) 
n=ni (I); 

t=sqrt ( (rA2) * (patm-patm/n) *sqrt (l-muA2) / (0 . 8*E) ) ;  %m  thickness  of 
sphere  based  on  E*  gives  t* 

ts=t*0.1128;  %from  LTAsphere  gives  tsmeared  value  vs  t* 
mb (I) = (rhos*ts*4*pi*r A2+ ( (4/3) *pi* (r- 

. 5*ts ) A3 ) *rhoair/n) / ( (4/3) *pi* ( (r+.5*ts) A3) *rhoair ) ;  %mass /buoyancy 
ratio 


end 

figure ( ) 

plot (1 . /ni,mb, 1 . /ni, ones (size(ni))*l,'k:') 

xlabel (' Fraction  of  Sea  Level  Pressure  Inside  Sphere') 

ylabel ( ' W/B ' ) 

%  title (' Effect  of  Decreasing  Internal  Pressure') 


o,  o, 
o  o 

%M= [pressure/ sea 
(m)  ] 

M=  [  1  1 

0.94211695 

0.887046632 

0.834542314 

0.784603997 

0.737132988 

0.692030595 

0.649198125 

0.608536886 

0.569948187 

0.533432026 

0.498791019 

0.466025167 

0.435035776 


level  pressure 


0 


0  . 

.  9529 

500 

0  . 

.  9075 

1000 

0  . 

.8638 

1500 

0  . 

.  8217 

2000 

0  . 

.7812 

2500 

0  . 

.7423 

3000 

0  . 

.7048 

3500 

0  . 

.  6689 

4000 

0  . 

.  6343 

4500 

0  . 

.  6012 

5000 

0  . 

.5694 

5500 

0  . 

.5389 

6000 

0  . 

.5096 

6500 

density/sea  level 


density,  elevation 
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0.405724155 

0.4817 

7000 

0.377991611 

0.4549 

7500 

0.351838145 

0.4292 

8000 

0.327165063 

0.4047 

8500 

0.303972366 

0.3813 

9000 

0.28206267 

0.3589 

9500 

0.261534666 

0.3376 

10000 

0.24219097 

0.3172 

10500 

0.224031582 

0.2978 

11000 

0.207056501 

0.2755 

11500 

0.191463114 

0.2546 

12000 

0.176955342 

0.2354 

12500 

0.163631878 

0.2176 

13000 

0.151295337 

0.2012 

13500 

0.139847027 

0 . 186 

14000 

0.129286948 

0.172 

14500 

0.119516408 

0 . 159 

15000 

0.110535406 

0 . 147 

15500 

0.102146558 

0 . 1359 

16000 

0.094468295 

0.1256 

16500 

0.087342709 

0 . 1162 

17000 

0.080750062 

0.1074 

17500 

0.074660745 

0.0993 

18000 

0.069035283 

0.09182 

18500 

0.063824328 

0.08489 

19000 

0.059018011 

0.0785 

19500 

0.054566987 

0.07258 

20000 

0.039940785 

0.05266 

22000 

0.029331359 

0.03832 

24000 

0.021593881 

0.02797 

26000 

0.01594868 

0.02047 

28000 

0.011813472 

0.01503 

30000] ; 

for  1=1 : length (M (:, 1 ) ) 

t3=sqrt  (  (r  A2 )  *  (patm*M  (1,1)  )  *sqrt  ( l-muA2 )  /  ( 0 . 8*E)  )  ; 
t3s=t3*0 . 1128;  %from  LTAsphere  gives  tsmeared  value  vs  t* 
t2  (I) =t3s ; 

mb2 (I) = (rhos*t2 (1) *4*pi*rA2) / ( (4/3) *pi* (rA3) *rhoair*M (1,2)); 
mb3 ( I ) = (rhos*t3s*4*pi*r A2 ) / ( (4/3) *pi* (rA3) *rhoair*M (1,2) ) ; 
end 

figure ( ) 

plot  (M  ( 1 :  4 , 3 ) *3.28084, mb3 (1:4),M(1:4,3) *3.28084,mb2 (1:4),M(1:4,3)*3.280 

84, ones (4, 1) *1,  'k:  '  ) 

xlabel (' Pressure  Altitude,  ft') 

ylabel ( ' W/B ' ) 

%  title (' Effect  of  Increasing  Designed  Deployment  Elevation') 
legend (' Designed  For  Specific  Altitude ',' Designed  for  Sea  Level 
Pressure ' ) 


A.12  Altitude  Effect  on  W/B  (LTAspherealt.m) 

%Capt  Metlen  LTAspherealt 
clc;  clear  all;  close  all; 


properties  of  air 
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da=1.225;  %kg/mA3,  density  of  air  at  sea  level 


%M= [pressure/sea  level  pressure,  density/sea  level  density,  elevation 
(m)  ] 

M=[same  as  Appendix  A.  11]; 
diameter=6 . 5  ; 

index=f ind (ro*2>diameter) ; index=index ( 1 )  ; 

diameter=ro (index) *2; 

mb= [0.9  0.8  0.5  0.25  0.09] ; 

for  1=1:5 

Mb ( : , I) =mb (I) *b (index) ./ (4/3*pi* (diameter/2) A3*M ( : , 2) *rho) ; 
end 

figure ( ) 

plot (M ( 1  rend,  3) *3 . 2 80 84 , Mb ( 1 : end, 1 )  ,M(1 rend,  3) *3 . 2 80 84 ,  Mb  ( 1 :  end,  2 )  ,M(1 
end,  3) *3. 28084, Mb (1: end, 3 ) , M ( 1 : end,  3) *3. 28084, Mb (1: end,  4 )  ,  M ( 1 :  end,  3 )  *3 
280 84 , Mb ( 1 : end, 5 ) )  %*3. 28084  to  convert  to  feet 
hold  on 

[x, y] =meshgrid (0:1000:65000,0.05:0.1:1.4) ; 
cv  =  [0  0 ]  ; 

[cont,h]  =  contourf (x, y, y-1, cv) ; 

%  Set  the  figure  Renderer  to  OpenGL,  which  supports  transparency 
set(gcf,  'Renderer',  'OpenGL'); 

%  Find  all  the  objects  that  are  children  of  the  contourgroup  that  have 
the  FaceAlpha  property 

al  =  f indob j (h,  '-property',  'FaceAlpha'); 
k  =  .3; 

%  Change  the  FaceAlpha  property,  which  will  change  the 
objects ' transparency 

set(al,  'FaceAlpha',  k, ' FaceColor ' , [ 0  0  l]);%red 

title (' Effect  of  Altitude  on  Weight/Buoyancy  Ratio') 
xlabel (' Density  Altitude,  ft') 
ylabel ( ' Weight/Buoyancy  Ratio') 

legend ( ' W/B=0 . 9 '  ,  'W/B=0 . 8 ' ,  ' W/B=0 .5',  ' W/B=0 . 25 ' ,  'W/B=0 . 09 ' ,  ' Location ' , 

NorthEastOutside ' ) 

axis ( [ 0 , 65000, 0.05, 1.2]) 

A.13  Optimization  of  Isogrid  (LTAisogrid  opt.m) 

%Capt  Metlen  LTAisogrid  opt 
clc;  clear  all;  close  all; 

%%  Optimization 

%  variables  x=[h  t  b  d] 

%  Set  your  upper  and  lower  bounds 

r=.33;  %  m,  radius  of  sphere 

scaling  =  [le-003  le-005  le-005  le-003] ; 

ub=[le-l  le-2  le-2  le-1 ]. /scaling; 

lb=[le-5  le-7  le-7  le-5] ./scaling; 


Set  your  initial  guess 
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xo—  [  1 . 4296e-003  3.2317e-005  1.4625e-005  1 . 2734e-003] ./scaling; 


%  Atmospheric  conditions 

rho_air=l . 225;  %  kg/mA3,  density  of  air 

P  atm=1500*101 . 325; %  Pa,  pressure  air 

%  Material  Properties 

E=895*10A9;  %  Pa,  skin  stiffness  UHM  carbon  fiber 
v=0.3;  %  Poisson's  ratio 

sc=1786e6;  %  Pa,  compressive  strength  epoxy 
Em=3.5*10A9;  %  Pa,  stiffness  epoxy 
thetai=0 :pi/50 :pi;  %  fiber  orientation  angles 
rho  skin=1522;  %  kg/mA3,  skin  density 

%  calculate  angle  and  Ex  of  each  ply  add  together  to  get  composite  Ex 
(Ec) 

for  1=1 : length (thetai) 
theta=thetai (I) ; 
thetal=theta+pi/2 ;  %90  degrees 
theta2=thetatpi/4 ;  %  45  degrees 
theta3=thetat3*pi/4 ;  %-45  degrees 

Ec ( I ) =0 . 6* (.25* (cos (theta) ) A4+.25* (cos (thetai ) ) A4+.25* (cos (theta2 ) ) A4+ . 

25* (cos ( theta3 ) ) A4 ) *E+ . 4  *Em; 

end 

E=E* . 6+Em* . 4 ; 

Eiso=min (Ec) ;  %min  fiber  stiffness  angle  for  quasi-isotropic  [ 9 0 / H — 
45/0] s 

%%  fmincon  Optimization 

%  This  section  runs  a  fmincon  optimization  routine 

const an ts= [r, rho_air, P^atm, E, Eiso, rho  skin, sc, v] ; 
options=optimset ( ' algorithm ', 'interior-point', 'display', 'iter'); 
[x,fval]  =  fmincon (@(x) 

LTAisogridcost (x, constants, scaling) , xo, [ ] , [ ] , [ ] , [ ] , lb, ub, @ (x) 
LTAisogrid_constraints (x, constants, scaling) , options) ; 

x=x . *scaling 

f=2*l . 05* . 33/ (sqrt (3) *x (1) )  %calculate  geometric  frequency 
f=round ( f ) 

%  Basic  information  on  geodesic  shape 
nV=10*fA2+2  %number  of  vertices 
nt=20*fA2  %number  of  triangles 
ne=30*fA2  %number  of  edges 


A.14  Optimization  of  Isogrid  Sub-Routine:  (LTAisogrid  constraints.m) 

%  Capt  Metlen 

function  [  c,ceq  ]  =  LTAisogrid_constraints(x, constants, scaling) 
x=x.*scaling; 
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%x=[h  t  b  d] 

%constants= [r, rho_air, P  atm, E, Eiso, rho^skin, sc, v] ; 

h=x ( 1 ) ; 
t=x ( 2 ) ; 
b=x ( 3 ) ; 
d=x ( 4 ) ; 

a=2*h/sqrt (3) ; 
t_bar=t+3*b*d/h; 

r=constants  (1) ; 

P_atm=constants  (3)  ; 

E=constants  (4)  ; 

Eiso=constants  (5)  ; 
sc=constants  (7) ; 
v=constants  (8)  ; 
bs= (E/Eiso) *b; 

alpha=bs*d/ (t*h) ; 
del=d/t; 

beta= (3*alpha* (1+del) A2+ (1+alpha) * ( l+alpha*del A2 ) ) A . 5; 
ts=t*beta/ (1+alpha) ; 

Es=Eiso* ( (1+alpha) A2) /beta; 

sigma=P  atm*r/ (2*t  bar) ;  %stress  in  skin  and  ribs 
sc_rib=0 . 456* (piA2) *E* ( (b/d) A2) / (l-vA2) ;  %rib  buckling  pressure 
sc_skin= ( 5/12 ) * (pi A2 ) *Eiso* ( ( t/a) A2 ) / ( l-vA2 ) ;  %skin  crippling  pressure 
P  crit=0 . 8 *Es* ( ( ts/r ) A2 ) /sqrt ( l-vA2 ) ;  %  global  sphere  buckling 

%constraints 
gl=sigma-sc; 
g2=sigma-sc_rib; 
g3=sigma-sc_skin; 
g4=P  atm-P  crit; 

c  =  [  g  1 ;  g  2 ;  g  3 ;  g  4  ]  ; 

ceq  =  [ ] ; 

end 


A.15  Optimization  of  Isogrid  Sub-Routine:  (LTAisogrid  cost.m) 

%  Capt  Metlen 

function  [  cost  ]  =  LTAisogrid_cost (  x, constants , scaling) 
x=x . *scaling; 

%x=[h  t  b  d] 

%constants= [r, rho_air , P, E, Eiso, rho_skin, sc, v] ; 

h=x ( 1 ) ; 
t=x ( 2 ) ; 
b=x ( 3 ) ; 
d=x ( 4 ) ; 

r=constants  (1)  ; 
rho_air=constants  (2)  ; 
rho  skin=constants  (  6 ) ; 

a  s=4*pi*rA2;  %  mA2,  surface  area  of  sphere 
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mA3,  volume  of  sphere 


v_s= (4/3) *pi*rA3;  % 
t_bar=t+3*b*d/h; 
cost  =  rho  skin*t  bar*a  s/ (rho_air*v  s);  %W/B  of  sphere 
end 
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virtual  elliptical  body  to  reduce  drag  and  maintain  structural  integrity  with  a  vacuum  was  also  investigated  and 
found  to  be  infeasible. 
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